diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 7ed7f090c0..2986024021 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -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 diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index 3bd29bb265..6dd104d5cd 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -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