From 33c52030a2a9d80afc2d6da0c3ccc4ed53d022c4 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 10 May 2016 05:32:04 -0400 Subject: [PATCH] Resolution of @paulromano comments --- src/mgxs_header.F90 | 194 ++++++++++++++-------------- src/scattdata_header.F90 | 269 ++++++++++++++++++++------------------- src/tally.F90 | 12 +- 3 files changed, 242 insertions(+), 233 deletions(-) diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 725aacffae..88c1b23e25 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -201,7 +201,7 @@ module mgxs_header ! the xsdata object node itself. !=============================================================================== - subroutine mgxs_init_file(this,node_xsdata,i_listing) + subroutine mgxs_init_file(this, node_xsdata, i_listing) class(Mgxs), intent(inout) :: this ! Working Object type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml integer, intent(in) :: i_listing ! Index in listings array @@ -236,7 +236,7 @@ module mgxs_header else if (temp_str == 'tabular') then this % scatt_type = ANGLE_TABULAR else - call fatal_error("Invalid Scatt Type Option!") + call fatal_error("Invalid scatt_type option!") end if else this % scatt_type = ANGLE_LEGENDRE @@ -259,8 +259,8 @@ module mgxs_header end subroutine mgxs_init_file - subroutine mgxsiso_init_file(this,node_xsdata,groups,get_kfiss,get_fiss, & - max_order,i_listing) + subroutine mgxsiso_init_file(this, node_xsdata, groups, get_kfiss, get_fiss, & + max_order, i_listing) class(MgxsIso), intent(inout) :: this ! Working Object type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml integer, intent(in) :: groups ! Number of Energy groups @@ -282,7 +282,7 @@ module mgxs_header integer :: legendre_mu_points, imu ! Call generic data gathering routine (will populate the metadata) - call mgxs_init_file(this,node_xsdata,i_listing) + call mgxs_init_file(this, node_xsdata, i_listing) ! Load the more specific data allocate(this % nu_fission(groups)) @@ -292,7 +292,7 @@ module mgxs_header ! Chi was provided, that means they are giving chi and nu-fission ! vectors ! Get chi - allocate(temp_arr(1 * groups)) + allocate(temp_arr(groups)) call get_node_array(node_xsdata, "chi", temp_arr) do gin = 1, groups do gout = 1, groups @@ -329,7 +329,7 @@ module mgxs_header end do ! Now pull out information needed for chi - this % chi = temp_2d + this % chi(:, :) = temp_2d ! Normalize chi so its CDF goes to 1 do gin = 1, groups this % chi(:, gin) = this % chi(:, gin) / sum(this % chi(:, gin)) @@ -375,14 +375,14 @@ module mgxs_header if (arr_len == groups * groups) then allocate(temp_arr(arr_len)) call get_node_array(node_xsdata, "multiplicity", temp_arr) - temp_mult = reshape(temp_arr, (/groups, groups/)) + temp_mult(:, :) = reshape(temp_arr, (/groups, groups/)) deallocate(temp_arr) else call fatal_error("Multiplicity length not same as number of groups& & squared!") end if else - temp_mult = ONE + temp_mult(:, :) = ONE end if ! Get scattering treatment information @@ -426,7 +426,7 @@ module mgxs_header if (check_for_node(node_xsdata, "order")) then call get_node_value(node_xsdata, "order", order) else - call fatal_error("Order Must Be Provided!") + call fatal_error("Order must be provided!") end if ! Before retrieving the data, store the dimensionality of the data in @@ -546,11 +546,11 @@ module mgxs_header if (check_for_node(node_xsdata, "total")) then call get_node_array(node_xsdata, "total", this % total) else - this % total = this % absorption + this % scatter % scattxs + this % total(:) = this % absorption(:) + this % scatter % scattxs(:) end if ! Deallocate temporaries for the next material - deallocate(input_scatt,scatt_coeffs,temp_mult) + deallocate(input_scatt, scatt_coeffs, temp_mult) ! Finally, check sigT to ensure it is not 0 since it is ! often divided by in the tally routines @@ -561,8 +561,8 @@ module mgxs_header end subroutine mgxsiso_init_file - subroutine mgxsang_init_file(this,node_xsdata,groups,get_kfiss,get_fiss, & - max_order,i_listing) + subroutine mgxsang_init_file(this, node_xsdata, groups, get_kfiss, get_fiss, & + max_order, i_listing) class(MgxsAngle), intent(inout) :: this ! Working Object type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml integer, intent(in) :: groups ! Number of Energy groups @@ -584,18 +584,18 @@ module mgxs_header integer :: legendre_mu_points, imu, ipol, iazi ! Call generic data gathering routine (will populate the metadata) - call mgxs_init_file(this,node_xsdata,i_listing) + call mgxs_init_file(this, node_xsdata, i_listing) if (check_for_node(node_xsdata, "num_polar")) then call get_node_value(node_xsdata, "num_polar", this % n_pol) else - call fatal_error("num_polar Must Be Provided!") + call fatal_error("num_polar must be provided!") end if if (check_for_node(node_xsdata, "num_azimuthal")) then call get_node_value(node_xsdata, "num_azimuthal", this % n_azi) else - call fatal_error("num_azimuthal Must Be Provided!") + call fatal_error("num_azimuthal must be provided!") end if ! Load angle data, if present (else equally spaced) @@ -663,8 +663,8 @@ module mgxs_header if (check_for_node(node_xsdata, "nu_fission")) then allocate(temp_arr(groups * this % n_azi * this % n_pol)) call get_node_array(node_xsdata, "nu_fission", temp_arr) - this % nu_fission = reshape(temp_arr,(/groups, this % n_azi, & - this % n_pol/)) + this % nu_fission(:, :, :) = reshape(temp_arr, (/groups, & + this % n_azi, this % n_pol/)) deallocate(temp_arr) else call fatal_error("If fissionable, must provide nu_fission!") @@ -677,8 +677,8 @@ module mgxs_header allocate(temp_arr(groups * groups * this % n_azi * this % n_pol)) call get_node_array(node_xsdata, "nu_fission", temp_arr) allocate(temp_4d(groups, groups, this % n_azi,this % n_pol)) - temp_4d = reshape(temp_arr, (/groups, groups, this % n_azi, & - this % n_pol/)) + temp_4d(:, :, :, :) = reshape(temp_arr, (/groups, groups, & + this % n_azi, this % n_pol/)) deallocate(temp_arr) else call fatal_error("If fissionable, must provide nu_fission!") @@ -716,8 +716,8 @@ module mgxs_header allocate(temp_arr(groups * this % n_azi * this % n_pol)) call get_node_array(node_xsdata, "fission", temp_arr) allocate(this % fission(groups, this % n_azi, this % n_pol)) - this % fission = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) + this % fission(:, :, :) = reshape(temp_arr, (/groups, this % n_azi, & + this % n_pol/)) deallocate(temp_arr) else call fatal_error("Fission data missing, required due to fission& @@ -729,8 +729,8 @@ module mgxs_header allocate(temp_arr(groups * this % n_azi * this % n_pol)) call get_node_array(node_xsdata, "kappa_fission", temp_arr) allocate(this % k_fission(groups, this % n_azi, this % n_pol)) - this % k_fission = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) + this % k_fission(:, :, :) = reshape(temp_arr, (/groups, & + this % n_azi, this % n_pol/)) deallocate(temp_arr) else call fatal_error("kappa_fission data missing, required due to & @@ -738,16 +738,16 @@ module mgxs_header end if end if else - this % nu_fission = ZERO - this % chi = ZERO + this % nu_fission(:, :, :) = ZERO + this % chi(:, :, :, :) = ZERO end if if (check_for_node(node_xsdata, "absorption")) then allocate(temp_arr(groups * this % n_azi * this % n_pol)) call get_node_array(node_xsdata, "absorption", temp_arr) allocate(this % absorption(groups, this % n_azi, this % n_pol)) - this % absorption = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) + this % absorption(:, :, :) = reshape(temp_arr, (/groups, this % n_azi, & + this % n_pol/)) deallocate(temp_arr) else call fatal_error("Must provide absorption!") @@ -760,15 +760,15 @@ module mgxs_header if (arr_len == groups * groups * this % n_azi * this % n_pol) then allocate(temp_arr(arr_len)) call get_node_array(node_xsdata, "multiplicity", temp_arr) - temp_mult = reshape(temp_arr, (/groups, groups, this % n_azi, & - this % n_pol/)) + temp_mult(:, :, :, :) = reshape(temp_arr, (/groups, groups, & + this % n_azi, this % n_pol/)) deallocate(temp_arr) else call fatal_error("Multiplicity length not same as number of groups& & squared!") end if else - temp_mult = ONE + temp_mult(:, :, :, :) = ONE end if ! Get scattering treatment information @@ -812,7 +812,7 @@ module mgxs_header if (check_for_node(node_xsdata, "order")) then call get_node_value(node_xsdata, "order", order) else - call fatal_error("Order Must Be Provided!") + call fatal_error("Order must be provided!") end if ! Before retrieving the data, store the dimensionality of the data in @@ -836,8 +836,8 @@ module mgxs_header allocate(temp_arr(groups * groups * order_dim * this % n_azi * & this % n_pol)) call get_node_array(node_xsdata, "scatter", temp_arr) - input_scatt = reshape(temp_arr, (/groups, groups, order_dim, & - this % n_azi, this % n_pol/)) + 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 @@ -951,8 +951,8 @@ module mgxs_header if (check_for_node(node_xsdata, "total")) then allocate(temp_arr(groups * this % n_azi * this % n_pol)) call get_node_array(node_xsdata, "total", temp_arr) - this % total = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) + this % total(:, :, :) = reshape(temp_arr, (/groups, this % n_azi, & + this % n_pol/)) deallocate(temp_arr) else do ipol = 1, this % n_pol @@ -1251,7 +1251,7 @@ module mgxs_header if (present(gout)) then xs = this % chi(gout, gin, iazi, ipol) else - ! Not sure youd want a 1 or a 0, but here you go! + ! Not sure you would want a 1 or a 0, but here you go! xs = sum(this % chi(:, gin, iazi, ipol)) end if case('scatter') @@ -1316,7 +1316,7 @@ module mgxs_header ! objects !=============================================================================== - subroutine mgxs_combine(this,mat,scatt_type,i_listing) + subroutine mgxs_combine(this, mat, scatt_type, i_listing) class(Mgxs), intent(inout) :: this ! The Mgxs to initialize type(Material), pointer, intent(in) :: mat ! base material integer, intent(in) :: scatt_type ! How is data presented @@ -1340,7 +1340,7 @@ module mgxs_header end subroutine mgxs_combine - subroutine mgxsiso_combine(this,mat,nuclides,groups,max_order,scatt_type, & + subroutine mgxsiso_combine(this, mat, nuclides, groups, max_order, scatt_type, & i_listing) class(MgxsIso), intent(inout) :: this ! The Mgxs to initialize type(Material), pointer, intent(in) :: mat ! base material @@ -1373,9 +1373,9 @@ module mgxs_header do i = 2, mat % n_nuclides select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsIso) - if (order /= size(nuc % scatter % dist(1) % data,dim=1)) & - call fatal_error("All Histogram Scattering Entries Must Be& - & Same Length!") + if (order /= size(nuc % scatter % dist(1) % data, dim=1)) & + call fatal_error("All histogram scattering entries must be& + & same length!") end select end do ! Ok, got our order, store the dimensionality @@ -1390,8 +1390,8 @@ module mgxs_header select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsIso) if (order /= size(nuc % scatter % dist(1) % data, dim=1)) & - call fatal_error("All Tabular Scattering Entries Must Be& - & Same Length!") + call fatal_error("All tabular scattering entries must be& + & same length!") end select end do ! Ok, got our order, store the dimensionality @@ -1406,7 +1406,7 @@ module mgxs_header do i = 1, mat % n_nuclides select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsIso) - if (size(nuc % scatter % dist(1) % data,dim=1) > mat_max_order) & + if (size(nuc % scatter % dist(1) % data, dim=1) > mat_max_order) & mat_max_order = size(nuc % scatter % dist(1) % data, dim=1) end select end do @@ -1423,25 +1423,25 @@ module mgxs_header ! Allocate and initialize data needed for macro_xs(i_mat) object allocate(this % total(groups)) - this % total = ZERO + this % total(:) = ZERO allocate(this % absorption(groups)) - this % absorption = ZERO + this % absorption(:) = ZERO allocate(this % fission(groups)) - this % fission = ZERO + this % fission(:) = ZERO allocate(this % k_fission(groups)) - this % k_fission = ZERO + this % k_fission(:) = ZERO allocate(this % nu_fission(groups)) - this % nu_fission = ZERO + this % nu_fission(:) = ZERO allocate(this % chi(groups,groups)) - this % chi = ZERO + this % chi(:, :) = ZERO allocate(temp_mult(groups,groups)) - temp_mult = ZERO + temp_mult(:, :) = ZERO allocate(mult_num(groups,groups)) - mult_num = ZERO + mult_num(:, :) = ZERO allocate(mult_denom(groups,groups)) - mult_denom = ZERO + mult_denom(:, :) = ZERO allocate(scatt_coeffs(order_dim,groups,groups)) - scatt_coeffs = ZERO + scatt_coeffs(:, :, :) = ZERO ! Add contribution from each nuclide in material do i = 1, mat % n_nuclides @@ -1452,19 +1452,19 @@ module mgxs_header select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsIso) ! Add contributions to total, absorption, and fission data (if necessary) - this % total = this % total + atom_density * nuc % total - this % absorption = this % absorption + & - atom_density * nuc % absorption + this % total(:) = this % total(:) + atom_density * nuc % total(:) + this % absorption(:) = this % absorption(:) + & + atom_density * nuc % absorption(:) if (nuc % fissionable) then - this % chi = this % chi + atom_density * nuc % chi - this % nu_fission = this % nu_fission + atom_density * & - nuc % nu_fission + this % chi(:, :) = this % chi(:, :) + atom_density * nuc % chi(:, :) + this % nu_fission(:) = this % nu_fission(:)+ atom_density * & + nuc % nu_fission(:) if (allocated(nuc % fission)) then - this % fission = this % fission + atom_density * nuc % fission + this % fission(:) = this % fission(:) + atom_density * nuc % fission(:) end if if (allocated(nuc % k_fission)) then - this % k_fission = this % k_fission + atom_density * & - nuc % k_fission + this % k_fission(:) = this % k_fission(:) + atom_density * & + nuc % k_fission(:) end if end if @@ -1498,7 +1498,7 @@ module mgxs_header nuc % scatter % get_matrix(min(nuc_order_dim, order_dim)) type is (MgxsAngle) - call fatal_error("Invalid Passing of MgxsAngle to MgxsIso Object") + call fatal_error("Invalid passing of MgxsAngle to MgxsIso object") end select end do @@ -1531,7 +1531,7 @@ module mgxs_header end subroutine mgxsiso_combine - subroutine mgxsang_combine(this,mat,nuclides,groups,max_order,scatt_type,& + subroutine mgxsang_combine(this, mat, nuclides, groups, max_order, scatt_type, & i_listing) class(MgxsAngle), intent(inout) :: this ! The Mgxs to initialize type(Material), pointer, intent(in) :: mat ! base material @@ -1551,7 +1551,7 @@ module mgxs_header real(8), allocatable :: mult_denom(:, :, :, :), scatt_coeffs(:, :, :, :, :) ! Set the meta-data - call mgxs_combine(this,mat,scatt_type,i_listing) + call mgxs_combine(this, mat, scatt_type, i_listing) ! Get the number of each polar and azi angles and make sure all the ! NuclideAngle types have the same number of these angles @@ -1564,12 +1564,12 @@ module mgxs_header n_pol = nuc % n_pol n_azi = nuc % n_azi allocate(this % polar(n_pol)) - this % polar = nuc % polar + this % polar(:) = nuc % polar(:) allocate(this % azimuthal(n_azi)) - this % azimuthal = nuc % azimuthal + this % azimuthal(:) = nuc % azimuthal(:) else if ((n_pol /= nuc % n_pol) .or. (n_azi /= nuc % n_azi)) then - call fatal_error("All Angular Data Must Be Same Length!") + call fatal_error("All angular data must be same length!") end if end if end select @@ -1589,8 +1589,8 @@ module mgxs_header select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsAngle) if (order /= size(nuc % scatter(1,1) % obj % dist(1) % data, dim=1)) & - call fatal_error("All Histogram Scattering Entries Must Be& - & Same Length!") + call fatal_error("All histogram scattering entries must be& + & same length!") end select end do ! Ok, got our order, store the dimensionality @@ -1609,9 +1609,9 @@ module mgxs_header do i = 2, mat % n_nuclides select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsAngle) - if (order /= size(nuc % scatter(1, 1) % obj % dist(1) % data,dim=1)) & - call fatal_error("All Tabular Scattering Entries Must Be& - & Same Length!") + if (order /= size(nuc % scatter(1, 1) % obj % dist(1) % data, dim=1)) & + call fatal_error("All tabular scattering entries must be& + & same length!") end select end do ! Ok, got our order, store the dimensionality @@ -1653,25 +1653,25 @@ module mgxs_header ! Allocate and initialize data within macro_xs(i_mat) object allocate(this % total(groups, n_azi, n_pol)) - this % total = ZERO + this % total(:, :, :) = ZERO allocate(this % absorption(groups, n_azi, n_pol)) - this % absorption = ZERO + this % absorption(:, :, :) = ZERO allocate(this % fission(groups, n_azi, n_pol)) - this % fission = ZERO + this % fission(:, :, :) = ZERO allocate(this % k_fission(groups, n_azi, n_pol)) - this % k_fission = ZERO + this % k_fission(:, :, :) = ZERO allocate(this % nu_fission(groups, n_azi, n_pol)) - this % nu_fission = ZERO + this % nu_fission(:, :, :) = ZERO allocate(this % chi(groups, groups, n_azi, n_pol)) - this % chi = ZERO + this % chi(:, :, :, :) = ZERO allocate(temp_mult(groups, groups, n_azi, n_pol)) - temp_mult = ZERO + temp_mult(:, :, :, :) = ZERO allocate(mult_num(groups, groups, n_azi, n_pol)) - mult_num = ZERO + mult_num(:, :, :, :) = ZERO allocate(mult_denom(groups, groups, n_azi, n_pol)) - mult_denom = ZERO + mult_denom(:, :, :, :) = ZERO allocate(scatt_coeffs(order_dim, groups, groups, n_azi, n_pol)) - scatt_coeffs = ZERO + scatt_coeffs(:, :, :, :, :) = ZERO ! Add contribution from each nuclide in material do i = 1, mat % n_nuclides @@ -1681,22 +1681,24 @@ module mgxs_header ! Perform our operations which depend upon the type select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsIso) - call fatal_error("Invalid Passing of MgxsIso to MgxsAngle Object") + call fatal_error("Invalid passing of MgxsIso to MgxsAngle object") type is (MgxsAngle) ! Add contributions to total, absorption, and fission data (if necessary) - this % total = this % total + atom_density * nuc % total - this % absorption = this % absorption + & - atom_density * nuc % absorption + this % total(:, :, :) = this % total(:, :, :) + & + atom_density * nuc % total(:, :, :) + this % absorption(:, :, :) = this % absorption(:, :, :) + & + atom_density * nuc % absorption(:, :, :) if (nuc % fissionable) then this % chi = this % chi + atom_density * nuc % chi - this % nu_fission = this % nu_fission + atom_density * & - nuc % nu_fission + this % nu_fission(:, :, :) = this % nu_fission(:, :, :) + & + atom_density * nuc % nu_fission(:, :, :) if (allocated(nuc % fission)) then - this % fission = this % fission + atom_density * nuc % fission + this % fission(:, :, :) = this % fission(:, :, :) + & + atom_density * nuc % fission(:, :, :) end if if (allocated(nuc % k_fission)) then - this % k_fission = this % k_fission + atom_density * & - nuc % k_fission + this % k_fission(:, :, :) = this % k_fission(:, :, :) + & + atom_density * nuc % k_fission(:, :, :) end if end if @@ -1730,7 +1732,7 @@ module mgxs_header end do ! Get the complete scattering matrix - nuc_order_dim = size(nuc % scatter(1,1) % obj % dist(1) % data,dim=1) + nuc_order_dim = size(nuc % scatter(1, 1) % obj % dist(1) % data, dim=1) do ipol = 1, n_pol do iazi = 1, n_azi scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) = & diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index f36fe6043e..12c11e2e4e 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -14,7 +14,7 @@ module scattdata_header !=============================================================================== type :: Jagged2D - real(8), allocatable :: data(:,:) + real(8), allocatable :: data(:, :) end type Jagged2D type :: Jagged1D @@ -27,12 +27,15 @@ module scattdata_header !=============================================================================== type, abstract :: ScattData - ! normalized p0 matrix on its own for sampling energy + ! The data attribute of the energy, mult, and dist arrays + ! are not necessarily 1-indexed as they instead will be allocated + ! from a minimum outgoing group to an outgoing minimum group. + ! Normalized p0 matrix on its own for sampling energy type(Jagged1D), allocatable :: energy(:) ! (Gin % data(Gout)) - ! nu-scatter multiplication (i.e. nu-scatt/scatt) + ! Nu-scatter multiplication (i.e. nu-scatt/scatt) type(Jagged1D), allocatable :: mult(:) ! (Gin % data(Gout)) ! Angular distribution - type(Jagged2D), allocatable :: dist(:) ! (Gin % data(Order/Nmu x Gout) + type(Jagged2D), allocatable :: dist(:) ! (Gin % data(Order/Nmu, Gout) integer, allocatable :: gmin(:) ! Minimum outgoing group integer, allocatable :: gmax(:) ! Maximum outgoing group real(8), allocatable :: scattxs(:) ! Isotropic Sigma_{s,g_{in}} @@ -47,9 +50,9 @@ module scattdata_header abstract interface subroutine scattdata_init_(this, mult, coeffs) import ScattData - class(ScattData), intent(inout) :: this ! Scattering Object to work with - real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix - real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use + 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_ pure function scattdata_calc_f_(this, gin, gout, mu) result(f) @@ -120,10 +123,10 @@ contains !=============================================================================== subroutine scattdata_init(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 + 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 integer :: groups, gmin, gmax, gin real(8) :: norm @@ -139,15 +142,15 @@ contains ! Also set energy values when doing it 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 + 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 + 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 + if (energy(gmax, gin) > ZERO) exit end do ! Treat the case of all zeros if (gmin > gmax) then @@ -156,10 +159,10 @@ contains ! 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) + 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)) + this % mult(gin) % data(gmin:gmax) = mult(gmin:gmax, gin) + allocate(this % dist(gin) % data(order, gmin:gmax)) this % dist(gin) % data = ZERO this % gmin(gin) = gmin this % gmax(gin) = gmax @@ -167,38 +170,38 @@ contains end subroutine scattdata_init subroutine scattdatalegendre_init(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 + 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 real(8) :: dmu, mu, f, norm integer :: imu, Nmu, gout, gin, groups, order - real(8), allocatable :: energy(:,:) - real(8), allocatable :: matrix(:,:,:) + real(8), allocatable :: energy(:, :) + real(8), allocatable :: matrix(:, :, :) - groups = size(coeffs,dim=3) - order = size(coeffs,dim=1) + groups = size(coeffs, dim=3) + order = size(coeffs, 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(order, groups, groups)) + matrix (:, :, :)= coeffs ! 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) + this % scattxs(:) = sum(matrix(1, :, :), dim=1) - allocate(energy(groups,groups)) - energy = ZERO + allocate(energy(groups, groups)) + energy(:, :) = ZERO ! 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 + norm = matrix(1, gout, gin) + energy(gout, gin) = norm if (norm /= ZERO) then - matrix(:,gout,gin) = matrix(:,gout,gin) / norm + matrix(:, gout, gin) = matrix(:, gout, gin) / norm end if end do end do @@ -209,16 +212,16 @@ contains ! 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) + this % dist(gin) % data(:, gout) = matrix(:, gout, gin) end do allocate(this % max_val(gin) % data(this % gmin(gin):this % gmax(gin))) - this % max_val(gin) % data = ZERO + this % max_val(gin) % data(:) = ZERO end do ! Step through the polynomial with fixed number of points to identify ! the maximal value. Nmu = 1001 - dmu = TWO / real(Nmu - 1,8) + dmu = TWO / real(Nmu - 1, 8) do gin = 1, groups do gout = this % gmin(gin), this % gmax(gin) do imu = 1, Nmu @@ -228,7 +231,7 @@ contains else if (imu == Nmu) then mu = ONE else - mu = -ONE + real(imu - 1,8) * dmu + mu = -ONE + real(imu - 1, 8) * dmu end if ! Calculate probability f = this % calc_f(gin,gout,mu) @@ -245,37 +248,37 @@ contains subroutine scattdatahistogram_init(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 + real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix + real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use integer :: imu, gin, gout, groups, order real(8) :: norm - real(8), allocatable :: energy(:,:) - real(8), allocatable :: matrix(:,:,:) + real(8), allocatable :: energy(:, :) + real(8), allocatable :: matrix(:, :, :) - groups = size(coeffs,dim=3) - order = size(coeffs,dim=1) + groups = size(coeffs, dim=3) + order = size(coeffs, 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(order, groups, groups)) + matrix(:, :, :) = coeffs ! 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) + this % scattxs(:) = sum(sum(matrix(:, :, :), dim=1), dim=1) - allocate(energy(groups,groups)) - energy = ZERO + allocate(energy(groups, groups)) + energy(:, :) = ZERO ! 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(:, gout, gin)) + energy(gout, gin) = norm if (norm /= ZERO) then - matrix(:,gout,gin) = matrix(:,gout,gin) / norm + matrix(:, gout, gin) = matrix(:, gout, gin) / norm end if end do end do @@ -283,10 +286,10 @@ contains call scattdata_init(this, order, energy, mult) allocate(this % mu(order)) - this % dmu = TWO / real(order,8) + this % dmu = TWO / real(order, 8) this % mu(1) = -ONE do imu = 2, order - this % mu(imu) = -ONE + real(imu - 1,8) * this % dmu + this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu end do ! Integrate this histogram so we can avoid rejection sampling while @@ -297,21 +300,23 @@ contains this % gmin(gin):this % gmax(gin))) do gout = this % gmin(gin), this % gmax(gin) ! Store the histogram - this % fmu(gin) % data(:,gout) = matrix(:,gout,gin) + this % fmu(gin) % data(:, gout) = matrix(:, gout, gin) ! Integrate the histogram - this % dist(gin) % data(1,gout) = this % dmu * matrix(1,gout,gin) + this % dist(gin) % data(1, gout) = & + this % dmu * matrix(1, gout, gin) do imu = 2, order - this % dist(gin) % data(imu,gout) = this % dmu * matrix(imu,gout,gin) + & - this % dist(gin) % data(imu - 1,gout) + this % dist(gin) % data(imu, gout) = & + this % dmu * matrix(imu, gout, gin) + & + this % dist(gin) % data(imu - 1, gout) end do ! Now make sure integral norms to zero - norm = this % dist(gin) % data(order,gout) + norm = this % dist(gin) % data(order, gout) if (norm > ZERO) then - this % fmu(gin) % data(:,gout) = & - this % fmu(gin) % data(:,gout) / norm - this % dist(gin) % data(:,gout) = & - this % dist(gin) % data(:,gout) / norm + this % fmu(gin) % data(:, gout) = & + this % fmu(gin) % data(:, gout) / norm + this % dist(gin) % data(:, gout) = & + this % dist(gin) % data(:, gout) / norm end if end do end do @@ -320,27 +325,27 @@ contains subroutine scattdatatabular_init(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 + real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix + real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use integer :: imu, gin, gout, groups, order real(8) :: norm - real(8), allocatable :: energy(:,:) - real(8), allocatable :: matrix(:,:,:) + real(8), allocatable :: energy(:, :) + real(8), allocatable :: matrix(:, :, :) - groups = size(coeffs,dim=3) - order = size(coeffs,dim=1) + groups = size(coeffs, dim=3) + order = size(coeffs, 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(order, groups, groups)) + matrix(:, :, :) = coeffs ! Build the angular distribution mu values allocate(this % mu(order)) - this % dmu = TWO / real(order - 1,8) + this % dmu = TWO / real(order - 1, 8) this % mu(1) = -ONE do imu = 2, order - 1 - this % mu(imu) = -ONE + real(imu - 1,8) * this % dmu + this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu end do this % mu(order) = ONE @@ -353,24 +358,24 @@ contains norm = ZERO do gout = 1, groups do imu = 2, order - norm = norm + HALF * this % dmu * (matrix(imu - 1,gout,gin) + & - matrix(imu,gout,gin)) + norm = norm + HALF * this % dmu * (matrix(imu - 1, gout, gin) + & + matrix(imu, gout, gin)) end do end do this % scattxs(gin) = norm end do - allocate(energy(groups,groups)) - energy = ZERO + allocate(energy(groups, groups)) + energy(:, :) = ZERO ! Build energy transfer probability matrix from data in matrix do gin = 1, groups do gout = 1, groups norm = ZERO do imu = 2, order norm = norm + HALF * this % dmu * & - (matrix(imu - 1,gout,gin) + matrix(imu,gout,gin)) + (matrix(imu - 1, gout, gin) + matrix(imu, gout, gin)) end do - energy(gout,gin) = norm + energy(gout, gin) = norm end do end do call scattdata_init(this, order, energy, mult) @@ -383,12 +388,12 @@ contains do gout = this % gmin(gin), this % 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(:, gout, gin) ! Force positivity do imu = 1, order - if (this % fmu(gin) % data(imu,gout) < ZERO) then - this % fmu(gin) % data(imu,gout) = ZERO + if (this % fmu(gin) % data(imu, gout) < ZERO) then + this % fmu(gin) % data(imu, gout) = ZERO end if end do @@ -397,27 +402,27 @@ contains norm = ZERO do imu = 2, order norm = norm + HALF * this % dmu * & - (this % fmu(gin) % data(imu - 1,gout) + & - this % fmu(gin) % data(imu,gout)) + (this % fmu(gin) % data(imu - 1, gout) + & + this % fmu(gin) % data(imu, gout)) end do if (norm > ZERO) then - this % fmu(gin) % data(:,gout) = & - this % fmu(gin) % data(:,gout) / norm + this % fmu(gin) % data(:, gout) = & + this % fmu(gin) % data(:, gout) / norm end if ! Now create CDF from fmu with trapezoidal rule - this % dist(gin) % data(1,gout) = ZERO + this % dist(gin) % data(1, gout) = ZERO do imu = 2, order - this % dist(gin) % data(imu,gout) = & - this % dist(gin) % data(imu - 1,gout) + & - HALF * this % dmu * (this % fmu(gin) % data(imu - 1,gout) + & - this % fmu(gin) % data(imu,gout)) + this % dist(gin) % data(imu, gout) = & + this % dist(gin) % data(imu - 1, gout) + & + HALF * this % dmu * (this % fmu(gin) % data(imu - 1, gout) + & + this % fmu(gin) % data(imu, gout)) end do ! Ensure we normalize to 1 still - norm = this % dist(gin) % data(order,gout) + norm = this % dist(gin) % data(order, gout) if (norm > ZERO) then - this % dist(gin) % data(:,gout) = & - this % dist(gin) % data(:,gout) / norm + this % dist(gin) % data(:, gout) = & + this % dist(gin) % data(:, gout) / norm end if end do end do @@ -438,7 +443,7 @@ contains if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then f = ZERO else - f = evaluate_legendre(this % dist(gin) % data(:,gout),mu) + f = evaluate_legendre(this % dist(gin) % data(:, gout), mu) end if end function scattdatalegendre_calc_f @@ -457,12 +462,12 @@ contains else ! Find mu bin if (mu == ONE) then - imu = size(this % fmu(gin) % data,dim=1) + imu = size(this % fmu(gin) % data, dim=1) else - imu = floor((mu + ONE)/ this % dmu + ONE) + imu = floor((mu + ONE) / this % dmu + ONE) end if - f = this % fmu(gin) % data(imu,gout) + f = this % fmu(gin) % data(imu, gout) end if end function scattdatahistogram_calc_f @@ -482,15 +487,15 @@ contains else ! Find mu bin if (mu == ONE) then - imu = size(this % fmu(gin) % data,dim=1) - 1 + imu = size(this % fmu(gin) % data, dim=1) - 1 else - imu = floor((mu + ONE)/ this % dmu + ONE) + imu = floor((mu + ONE) / this % dmu + ONE) end if ! Now interpolate to find f(mu) r = (mu - this % mu(imu)) / (this % mu(imu + 1) - this % mu(imu)) - f = (ONE - r) * this % fmu(gin) % data(imu,gout) + & - r * this % fmu(gin) % data(imu + 1,gout) + f = (ONE - r) * this % fmu(gin) % data(imu, gout) + & + r * this % fmu(gin) % data(imu + 1, gout) end if end function scattdatatabular_calc_f @@ -529,7 +534,7 @@ contains samples = 0 do mu = TWO * prn() - ONE - f = this % calc_f(gin,gout,mu) + f = this % calc_f(gin, gout, mu) if (f > ZERO) then u = prn() * M if (u <= f) then @@ -567,11 +572,11 @@ contains end do xi = prn() - if (xi < this % dist(gin) % data(1,gout)) then + if (xi < this % dist(gin) % data(1, gout)) then imu = 1 else - imu = binary_search(this % dist(gin) % data(:,gout), & - size(this % dist(gin) % data(:,gout)), xi) + imu = binary_search(this % dist(gin) % data(:, gout), & + size(this % dist(gin) % data(:, gout)), xi) end if ! Randomly select a mu in this bin. @@ -604,12 +609,12 @@ contains end do ! determine outgoing cosine bin - NP = size(this % dist(gin) % data(:,gout)) + NP = size(this % dist(gin) % data(:, gout)) xi = prn() - c_k = this % dist(gin) % data(1,gout) + c_k = this % dist(gin) % data(1, gout) do k = 1, NP - 1 - c_k1 = this % dist(gin) % data(k + 1,gout) + c_k1 = this % dist(gin) % data(k + 1, gout) if (xi < c_k1) exit c_k = c_k1 end do @@ -617,18 +622,19 @@ contains ! check to make sure k is <= NP - 1 k = min(k, NP - 1) - p0 = this % fmu(gin) % data(k,gout) + p0 = this % fmu(gin) % data(k, gout) mu0 = this % mu(k) ! Linear-linear interpolation to find mu value w/in bin. - p1 = this % fmu(gin) % data(k + 1,gout) + p1 = this % fmu(gin) % data(k + 1, gout) mu1 = this % mu(k + 1) - frac = (p1 - p0)/(mu1 - mu0) + frac = (p1 - p0) / (mu1 - mu0) if (frac == ZERO) then - mu = mu0 + (xi - c_k)/p0 + mu = mu0 + (xi - c_k) / p0 else - mu = mu0 + (sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac + mu = mu0 + & + (sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac end if if (mu <= -ONE) then @@ -649,23 +655,24 @@ contains 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 + real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built integer :: order, groups, gin, gout groups = size(this % energy) - order = min(req_order,size(this % dist(1) % data(:,1))) + ! Set gin and gout for getting the order + order = min(req_order, size(this % dist(1) % data, dim=1)) - allocate(matrix(order,groups,groups)) + allocate(matrix(order, groups, 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 + matrix(:, :, :) = ZERO do gin = 1, groups do gout = this % gmin(gin), this % gmax(gin) - matrix(:,gout,gin) = this % scattxs(gin) * & + matrix(:, gout, gin) = this % scattxs(gin) * & this % energy(gin) % data(gout) * & - this % dist(gin) % data(1:order,gout) + this % dist(gin) % data(1:order, gout) end do end do end function scattdata_get_matrix @@ -673,23 +680,23 @@ contains 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 + real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built integer :: order, groups, gin, gout groups = size(this % energy) - order = min(req_order,size(this % dist(1) % data(:,1))) + order = min(req_order, size(this % dist(1) % data, dim=1)) - allocate(matrix(order,groups,groups)) + allocate(matrix(order, groups, 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 + matrix(:, :, :) = ZERO do gin = 1, groups do gout = this % gmin(gin), this % gmax(gin) - matrix(:,gout,gin) = this % scattxs(gin) * & + matrix(:, gout, gin) = this % scattxs(gin) * & this % energy(gin) % data(gout) * & - this % fmu(gin) % data(1:order,gout) + this % fmu(gin) % data(1:order, gout) end do end do end function scattdatahistogram_get_matrix @@ -697,23 +704,23 @@ contains 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 + real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built integer :: order, groups, gin, gout groups = size(this % energy) - order = min(req_order,size(this % dist(1) % data(:,1))) + order = min(req_order, size(this % dist(1) % data, dim=1)) - allocate(matrix(order,groups,groups)) + allocate(matrix(order, groups, 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 + matrix(:, :, :) = ZERO do gin = 1, groups do gout = this % gmin(gin), this % gmax(gin) - matrix(:,gout,gin) = this % scattxs(gin) * & + matrix(:, gout, gin) = this % scattxs(gin) * & this % energy(gin) % data(gout) * & - this % fmu(gin) % data(1:order,gout) + this % fmu(gin) % data(1:order, gout) end do end do end function scattdatatabular_get_matrix diff --git a/src/tally.F90 b/src/tally.F90 index 86e108c3dc..c4eaf30c8c 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -961,9 +961,9 @@ contains if (i_nuclide > 0) then score = score * atom_density * & nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, & - UVW=p_uvw,MU=p % mu) / & + UVW=p_uvw, MU=p % mu) / & matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, & - UVW=p_uvw,MU=p % mu) + UVW=p_uvw, MU=p % mu) end if else @@ -1080,11 +1080,11 @@ contains end if if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & + nucxs % get_xs('fission', p_g, UVW=p_uvw) / & matxs % get_xs('absorption', p_g, UVW=p_uvw) else score = score * & - matxs % get_xs('fission', p_g, UVW=p_uvw) / & + matxs % get_xs('fission', p_g, UVW=p_uvw) / & matxs % get_xs('absorption', p_g, UVW=p_uvw) end if else @@ -1170,11 +1170,11 @@ contains if (i_nuclide > 0) then score = score * atom_density * & nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + matxs % get_xs('absorption', p_g, UVW=p_uvw) else score = score * & matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + matxs % get_xs('absorption', p_g, UVW=p_uvw) end if else if (i_nuclide > 0) then