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

Updated fortran side of the sparse scattering matrices, almost

Browse source
This commit is contained in:
Adam Nelson 2016-09-07 18:39:10 -04:00
parent 1ef9acb8bd
commit 0cc5f06114
4 changed files with 193 additions and 132 deletions
Download patch file Download diff file Expand all files Collapse all files

42
docs/source/io_formats/mgxs_library.rst
View file

@ -118,20 +118,27 @@ Temperature-dependent data, provided for temperature <TTT>K.
Data specific to neutron scattering for the temperature <TTT>K
:Datasets: - **g_out bounds** (*int[2]* or *int[][][2]) --
Minimum (most energetic) and maximum (most thermal) outgoing groups
with non-zero values of the scattering matrix. These group numbers
use the standard ordering where the fastest neutron energy group
is group 1 while the most thermal neutron energy group is group G.
:Datasets: - **g_min** (*int[]* or *int[][][]) --
Minimum (most energetic) outgoing groups with non-zero values of
the scattering matrix. These group numbers use the standard
ordering where the fastest neutron energy group is group 1 while
the slowest neutron energy group is group G.
The dimensionality of `g_out bounds` is:
`g_out bounds][g_in][g_out]`, or
`g_out bounds[num-polar][num-azimuthal][g_in][g_out]`.
`g_min[g_in]`, or `g_min[num-polar][num-azimuthal][g_in]`.
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle".
- **scatter matrix** (*double[][]* or *double[][][][]*) --
Flattened representation of the scattering moment matrices where
the [g_in] and [g_out] indices are flattened. The pre-flattened
array is shaped as follows (in row-major format):
- **g_max** (*int[]* or *int[][][]) --
Maximum (least energetic) outgoing groups with non-zero values of
the scattering matrix. These group numbers use the standard
ordering where the fastest neutron energy group is group 1 while
the slowest neutron energy group is group G.
The dimensionality of `g_out bounds` is:
`g_max[g_in]`, or `g_max[num-polar][num-azimuthal][g_in]`.
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle".
- **scatter matrix** (*double[]*) -- Flattened representation of the
scattering moment matrices. The pre-flattened array is shaped as
follows (in row-major format):
`scatter matrix[order(+1)][g_in][g_out]`, or
`scatter matrix[num-polar][num-azimuthal][order(+1)][g_in][g_out]`
The former is used when `representation` is "isotropic", and the
@ -141,13 +148,12 @@ Data specific to neutron scattering for the temperature <TTT>K
Finally, the g_out dimension has a dimensionality of
`g_out bounds`[0] to `g_out bounds`[1].
- **multiplicity matrix** (*double[]*) -- Flattened representation of
the scattering moment matrices where the [g_in] and [g_out] indices
are flattened. This dataset provides the code with a scaling factor
to account for neutrons being produced in (n,xn) reactions. This
is assumed isotropic and therefore is not repeated for every
Legendre moment or histogram/tabular bin. This dataset is optional,
if it is not provided no multiplication (i.e., values of 1.0) will
be assumed.
the scattering moment matrices. This dataset provides the code with
a scaling factor to account for neutrons being produced in (n,xn)
reactions. This is assumed isotropic and therefore is not repeated
for every Legendre moment or histogram/tabular bin. This dataset is
optional, if it is not provided no multiplication (i.e., values of
1.0) will be assumed.
The pre-flattened array is shaped as follows (in row-major format):
`multiplicity matrix[g_in][g_out]`, or
`multiplicity matrix[num-polar][num-azimuthal][g_in][g_out]`

16
openmc/mgxs_library.py
View file

@ -1225,7 +1225,9 @@ class XSdata(object):
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :] += 1
scatt_grp.create_dataset("g_out bounds", data=g_out_bounds,
scatt_grp.create_dataset("g_min", data=g_out_bounds[:, 0],
compression=compression)
scatt_grp.create_dataset("g_max", data=g_out_bounds[:, 1],
compression=compression)
else:
@ -1240,7 +1242,7 @@ class XSdata(object):
g_out_bounds[p, a, g_in, 0] = nz[0]
g_out_bounds[p, a, g_in, 1] = nz[-1]
# Now create the flattened scatter matrix array
flat_scatt = [[[] for a in range(Na)] for p in range(Np)]
flat_scatt = []
for p in range(Np):
for a in range(Na):
for l in range(self._scatter_matrix[i].shape[2]):
@ -1250,16 +1252,16 @@ class XSdata(object):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1]
+ 1):
flat_scatt[p][a].append(matrix[g_out])
flat_scatt.append(matrix[g_out])
# And write it.
scatt_grp = xsgrp.create_group('scatter data')
scatt_grp.create_dataset("scatter matrix",
data=np.array(flat_scatt).flatten(),
data=np.array(flat_scatt),
compression=compression)
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
flat_mult = [[[] for a in range(Na)] for p in range(Np)]
flat_mult = []
for p in range(Np):
for a in range(Na):
for l in range(self._scatter_matrix[i].shape[2]):
@ -1268,9 +1270,9 @@ class XSdata(object):
self._multiplicity_matrix[i][p, a, g_in, :]
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
flat_mult[p][a].append(matrix[g_out])
flat_mult.append(matrix[g_out])
scatt_grp.create_dataset("multiplicity matrix",
data=np.array(flat_mult).flatten(),
data=np.array(flat_mult),
compression=compression)
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.

217
src/mgxs_header.F90
View file

@ -381,16 +381,17 @@ module mgxs_header
! in that conversion
character(MAX_LINE_LEN) :: temp_str
integer(HID_T) :: xsdata_grp
integer(HID_T) :: xsdata_grp, scatt_grp
real(8), allocatable :: temp_arr(:), temp_2d(:, :)
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
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(:)
! Call generic data gathering routine (will populate the metadata)
call mgxs_from_hdf5(this, xs_id, temperature, method, tolerance, &
@ -488,84 +489,134 @@ module mgxs_header
end if
! Get scattering data
! The input is gathered in the more user-friendly facing format of
! Gout x Gin x Order. 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).
allocate(input_scatt(groups, groups, order_dim))
if (check_dataset(xsdata_grp, "scatter matrix")) then
call read_dataset(input_scatt, xsdata_grp, "scatter matrix")
! 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))
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))
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) = &
evaluate_legendre(temp_scatt(gout, gin, :),mu)
! Ensure positivity of distribution
if (scatt_coeffs(imu, gout, gin) < ZERO) &
scatt_coeffs(imu, gout, gin) = ZERO
! And accrue the integral
if (imu > 1) then
norm = norm + HALF * dmu * &
(scatt_coeffs(imu - 1, gout, gin) + &
scatt_coeffs(imu, gout, gin))
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(:, gout, gin) = scatt_coeffs(:, gout, gin) * &
temp_scatt(gout, gin, 1) / norm
end if
end do ! gout
end do ! gin
else
! Sticking with current representation, carry forward but change
! the array ordering
allocate(scatt_coeffs(order_dim, groups, groups))
do gin = 1, groups
do gout = 1, groups
do l = 1, order_dim
scatt_coeffs(l, gout, gin) = temp_scatt(gout, gin, l)
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(g_min(groups))
call read_dataset(g_min, 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(g_max(groups))
call read_dataset(g_max, 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 gin = 1, groups
length = length + order_dim * (g_max(gin) - gmin(gin) + 1)
end do
! Allocate flattened array
allocate(temp_arr(length))
jf (.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))
index = 1
do gin = 1, groups
allocate(input_scatt(gin) % data(order_dim, gmin(gin):gmax(gin)))
do l = 1, order_dim
do gout = gmin(gin), gmax(gin)
input_scatt(gin) % data(l, gout) = temp_arr(index)
index = index + 1
end do
end do
end do
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))
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 gin = 1, groups
allocate(scatt_coeffs(gin) % data(order_dim, groups))
do gout = gmin(gin), gmax(gin)
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) % data(imu, gout) = &
evaluate_legendre(input_scatt(gin) % data(:, gout), mu)
! Ensure positivity of distribution
if (scatt_coeffs(gin) % data(imu, gout) < ZERO) &
scatt_coeffs(gin) % data(imu, gout) = ZERO
! And accrue the integral
if (imu > 1) then
norm = norm + HALF * dmu * &
(scatt_coeffs(gin) % data(imu - 1, gout) + &
scatt_coeffs(gin) % 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) % data(i:, gout) = &
scatt_coeffs(gin) % data(:, gout) * &
input_scatt(gin) % data(1, gout)
end if
end do ! gout
end do ! gin
else
! Sticking with current representation, carry forward but change
! the array ordering
do gin = 1, groups
allocate(scatt_coeffs(gin) % data(order_dim, groups))
scatt_coeffs(gin) % data(:, :) = input_scatt(gin) % data(:, :)
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 gin = 1, groups
length = length + (g_max(gin) - gmin(gin) + 1)
end do
! Allocate flattened array
allocate(temp_arr(length))
jf (.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))
index = 1
do gin = 1, groups
allocate(temp_mult(gin) % data(gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
temp_mult(gin) % data(gout) = temp_arr(index)
index = index + 1
end do
end do
deallocate(temp_arr)
! Allocate and initialize our ScattData Object.
if (this % scatter_type == ANGLE_HISTOGRAM) then
@ -577,7 +628,7 @@ module mgxs_header
end if
! Initialize the ScattData Object
call xs % scatter % init(temp_mult, scatt_coeffs)
call xs % scatter % init(gmin, gmax, temp_mult, scatt_coeffs)
! Check sigA to ensure it is not 0 since it is
! often divided by in the tally routines
@ -1168,7 +1219,7 @@ module mgxs_header
end do
! Initialize the ScattData Object
call this % xs(t) % scatter % init(temp_mult, scatt_coeffs)
call this % xs(t) % scatter % init_from_dense(temp_mult, scatt_coeffs)
! Now normalize chi
if (mat % fissionable) then
@ -1391,7 +1442,7 @@ module mgxs_header
! Initialize the ScattData Object
do ipol = 1, n_pol
do iazi = 1, n_azi
call this % xs(t) % scatter(iazi, ipol) % obj % init(&
call this % xs(t) % scatter(iazi, ipol) % obj % init_from_dense(&
temp_mult(:, :, iazi, ipol), &
scatt_coeffs(:, :, :, iazi, ipol))
end do

50
src/scattdata_header.F90
View file

@ -42,18 +42,20 @@ module scattdata_header
contains
procedure(scattdata_init_), deferred :: init ! Initializes ScattData
! Initializes ScattData from a dense matrix
procedure(scattdata_init_dense_), deferred :: init_from_dense
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_(this, mult, coeffs)
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_
end subroutine scattdata_init_dense_
pure function scattdata_calc_f_(this, gin, gout, mu) result(f)
import ScattData
@ -79,9 +81,9 @@ module scattdata_header
! Maximal value for rejection sampling from rectangle
type(Jagged1D), allocatable :: max_val(:) ! (Gin % data(Gout))
contains
procedure :: init => scattdatalegendre_init
procedure :: calc_f => scattdatalegendre_calc_f
procedure :: sample => scattdatalegendre_sample
procedure :: init_from_dense => scattdatalegendre_init_from_dense
procedure :: calc_f => scattdatalegendre_calc_f
procedure :: sample => scattdatalegendre_sample
end type ScattDataLegendre
type, extends(ScattData) :: ScattDataHistogram
@ -90,10 +92,10 @@ module scattdata_header
! Histogram of f(mu) (dist has CDF)
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
contains
procedure :: init => scattdatahistogram_init
procedure :: calc_f => scattdatahistogram_calc_f
procedure :: sample => scattdatahistogram_sample
procedure :: get_matrix => scattdatahistogram_get_matrix
procedure :: init_from_dense => scattdatahistogram_init_from_dense
procedure :: calc_f => scattdatahistogram_calc_f
procedure :: sample => scattdatahistogram_sample
procedure :: get_matrix => scattdatahistogram_get_matrix
end type ScattDataHistogram
type, extends(ScattData) :: ScattDataTabular
@ -102,10 +104,10 @@ module scattdata_header
! PDF of f(mu) (dist has CDF)
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
contains
procedure :: init => scattdatatabular_init
procedure :: calc_f => scattdatatabular_calc_f
procedure :: sample => scattdatatabular_sample
procedure :: get_matrix => scattdatatabular_get_matrix
procedure :: init_from_dense => scattdatatabular_init_from_dense
procedure :: calc_f => scattdatatabular_calc_f
procedure :: sample => scattdatatabular_sample
procedure :: get_matrix => scattdatatabular_get_matrix
end type ScattDataTabular
!===============================================================================
@ -122,7 +124,7 @@ contains
! SCATTDATA*_INIT builds the scattdata object
!===============================================================================
subroutine scattdata_init(this, order, energy, mult)
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
@ -167,9 +169,9 @@ contains
this % gmin(gin) = gmin
this % gmax(gin) = gmax
end do
end subroutine scattdata_init
end subroutine scattdata_init_from_dense
subroutine scattdatalegendre_init(this, mult, coeffs)
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
@ -206,7 +208,7 @@ contains
end do
end do
call scattdata_init(this, order, energy, mult)
call scattdata_init_from_dense(this, order, energy, mult)
allocate(this % max_val(groups))
! Set dist values from matrix and initialize max_val
@ -244,9 +246,9 @@ contains
this % max_val(gin) % data(gout) * 1.1_8
end do
end do
end subroutine scattdatalegendre_init
end subroutine scattdatalegendre_init_from_dense
subroutine scattdatahistogram_init(this, mult, coeffs)
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
@ -283,7 +285,7 @@ contains
end do
end do
call scattdata_init(this, order, energy, mult)
call scattdata_init_from_dense(this, order, energy, mult)
allocate(this % mu(order))
this % dmu = TWO / real(order, 8)
@ -321,9 +323,9 @@ contains
end do
end do
end subroutine scattdatahistogram_init
end subroutine scattdatahistogram_init_from_dense
subroutine scattdatatabular_init(this, mult, coeffs)
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
@ -378,7 +380,7 @@ contains
energy(gout, gin) = norm
end do
end do
call scattdata_init(this, order, energy, mult)
call scattdata_init_from_dense(this, order, energy, mult)
! Calculate f(mu) and integrate it so we can avoid rejection sampling
allocate(this % fmu(groups))
@ -415,7 +417,7 @@ contains
end if
end do
end do
end subroutine scattdatatabular_init
end subroutine scattdatatabular_init_from_dense
!===============================================================================
! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)