From 0cc5f06114aca33af769bc6318a74ca9b64f2c69 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 7 Sep 2016 18:39:10 -0400 Subject: [PATCH] Updated fortran side of the sparse scattering matrices, almost --- docs/source/io_formats/mgxs_library.rst | 42 +++-- openmc/mgxs_library.py | 16 +- src/mgxs_header.F90 | 217 +++++++++++++++--------- src/scattdata_header.F90 | 50 +++--- 4 files changed, 193 insertions(+), 132 deletions(-) diff --git a/docs/source/io_formats/mgxs_library.rst b/docs/source/io_formats/mgxs_library.rst index 4f6233aafc..0a02431f6c 100644 --- a/docs/source/io_formats/mgxs_library.rst +++ b/docs/source/io_formats/mgxs_library.rst @@ -118,20 +118,27 @@ Temperature-dependent data, provided for temperature K. Data specific to neutron scattering for the temperature 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 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]` diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 687e2d693a..eaf66d1200 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -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. diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 73a11476c5..7ed7f090c0 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -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 diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index d4643a0721..3bd29bb265 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -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)