diff --git a/src/ace.F90 b/src/ace.F90 index caa7c3aed0..4a9f287059 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -55,7 +55,7 @@ contains character(12) :: name ! name of isotope, e.g. 92235.03c character(12) :: alias ! alias of nuclide, e.g. U-235.03c type(Material), pointer :: mat - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc type(SAlphaBeta), pointer :: sab type(SetChar) :: already_read @@ -266,7 +266,7 @@ contains character(10) :: mat ! material identifier character(70) :: comment ! comment for ACE table character(MAX_FILE_LEN) :: filename ! path to ACE cross section library - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc type(SAlphaBeta), pointer :: sab type(XsListing), pointer :: listing @@ -423,7 +423,7 @@ contains !=============================================================================== subroutine read_esz(nuc, data_0K) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc logical, intent(in) :: data_0K ! are we reading 0K data? integer :: NE ! number of energy points for total and elastic cross sections @@ -511,7 +511,7 @@ contains !=============================================================================== subroutine read_nu_data(nuc) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i, j ! loop index integer :: idx ! index in XSS @@ -726,7 +726,7 @@ contains !=============================================================================== subroutine read_reactions(nuc) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i ! loop indices integer :: i_fission ! index in nuc % index_fission @@ -902,7 +902,7 @@ contains !=============================================================================== subroutine read_angular_dist(nuc) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc integer :: LOCB ! location of angular distribution for given MT integer :: NE ! number of incoming energies @@ -1006,7 +1006,7 @@ contains !=============================================================================== subroutine read_energy_dist(nuc) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i ! loop index integer :: n @@ -1395,7 +1395,7 @@ contains !=============================================================================== subroutine read_unr_res(nuc) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc integer :: JXS23 ! location of URR data integer :: lc ! locator @@ -1482,7 +1482,7 @@ contains !=============================================================================== subroutine generate_nu_fission(nuc) - type(NuclideCE), intent(inout) :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i ! index on nuclide energy grid diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b2813bef0c..fd8318c3b9 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -146,7 +146,7 @@ contains integer :: i_low ! lower logarithmic mapping index integer :: i_high ! upper logarithmic mapping index real(8) :: f ! interp factor on nuclide energy grid - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc type(Material), pointer :: mat ! Set pointer to nuclide and material @@ -513,9 +513,9 @@ contains !=============================================================================== pure function elastic_xs_0K(E, nuc) result(xs_out) - real(8), intent(in) :: E ! trial energy - type(NuclideCE), intent(in) :: nuc ! target nuclide at temperature - real(8) :: xs_out ! 0K xs at trial energy + real(8), intent(in) :: E ! trial energy + type(Nuclide), intent(in) :: nuc ! target nuclide at temperature + real(8) :: xs_out ! 0K xs at trial energy integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 248462f70c..66419f83cc 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -27,7 +27,7 @@ contains integer :: i ! index in nuclides array integer :: j ! index in materials array type(ListReal) :: list - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc type(Material), pointer :: mat call write_message("Creating unionized energy grid...", 5) @@ -70,7 +70,7 @@ contains real(8) :: E_max ! Maximum energy in MeV real(8) :: E_min ! Minimum energy in MeV real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc ! Set minimum/maximum energies E_max = energy_max_neutron @@ -179,7 +179,7 @@ contains integer :: index_e ! index on union energy grid real(8) :: union_energy ! energy on union grid real(8) :: energy ! energy on nuclide grid - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc type(Material), pointer :: mat do k = 1, n_materials diff --git a/src/global.F90 b/src/global.F90 index 6befa828ff..bb7bf50efc 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -5,9 +5,9 @@ module global use constants use dict_header, only: DictCharInt, DictIntInt use geometry_header, only: Cell, Universe, Lattice, LatticeContainer - use macroxs_header, only: MacroXSContainer use material_header, only: Material use mesh_header, only: RegularMesh + use mgxs_header, only: Mgxs, MgxsContainer use nuclide_header use plot_header, only: ObjectPlot use sab_header, only: SAlphaBeta @@ -86,7 +86,7 @@ module global ! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES ! Cross section arrays - type(NuclideCE), allocatable, target :: nuclides(:) ! Nuclide cross-sections + type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables @@ -112,10 +112,10 @@ module global ! MULTI-GROUP CROSS SECTION RELATED VARIABLES ! Cross section arrays - type(NuclideMGContainer), allocatable, target :: nuclides_MG(:) + type(MgxsContainer), allocatable, target :: nuclides_MG(:) ! Cross section caches - type(MacroXSContainer), target, allocatable :: macro_xs(:) + type(MgxsContainer), target, allocatable :: macro_xs(:) ! Number of energy groups integer :: energy_groups diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b69235690b..5a4bc8429d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -173,8 +173,12 @@ contains if (check_for_node(doc, "max_order")) then call get_node_value(doc, "max_order", max_order) else - ! Set to default of largest int, which means to use whatever is contained in library - max_order = huge(0) + ! Set to default of largest int - 1, which means to use whatever is + ! contained in library. + ! This is largest int - 1 because for legendre scattering, a value of + ! 1 is added to the order; adding 1 to huge(0) gets you the largest + ! negative integer, which is not what we want. + max_order = huge(0) - 1 end if else max_order = 0 @@ -2857,11 +2861,15 @@ contains allocate(t % filters(j) % real_bins(n_words)) call get_node_array(node_filt, "bins", t % filters(j) % real_bins) + ! We can save tallying time if we know that the tally bins + ! match the energy group structure. In that case, the matching bin + ! index is simply the group (after flipping for the different + ! ordering of the library and tallying systems). if (.not. run_CE) then - if (n_words /= energy_groups + 1) then - t % energy_matches_groups = .false. - else if (all(t % filters(j) % real_bins == energy_bins)) then - t % energy_matches_groups = .false. + if (n_words == energy_groups + 1) then + if (all(t % filters(j) % real_bins == & + energy_bins(energy_groups + 1:1:-1))) & + t % energy_matches_groups = .true. end if end if @@ -2876,11 +2884,15 @@ contains allocate(t % filters(j) % real_bins(n_words)) call get_node_array(node_filt, "bins", t % filters(j) % real_bins) + ! We can save tallying time if we know that the tally bins + ! match the energy group structure. In that case, the matching bin + ! index is simply the group (after flipping for the different + ! ordering of the library and tallying systems). if (.not. run_CE) then - if (n_words /= energy_groups + 1) then - t % energy_matches_groups = .false. - else if (all(t % filters(j) % real_bins == energy_bins)) then - t % energy_matches_groups = .false. + if (n_words == energy_groups + 1) then + if (all(t % filters(j) % real_bins == & + energy_bins(energy_groups + 1:1:-1))) & + t % energyout_matches_groups = .true. end if end if @@ -3338,7 +3350,13 @@ contains case ('nu-scatter') t % score_bins(j) = SCORE_NU_SCATTER - t % estimator = ESTIMATOR_ANALOG + + ! Set tally estimator to analog for CE mode + ! (MG mode has all data available without a collision being + ! necessary) + if (run_CE) then + t % estimator = ESTIMATOR_ANALOG + end if case ('scatter-n') t % score_bins(j) = SCORE_SCATTER_N @@ -3387,10 +3405,14 @@ contains call fatal_error("Diffusion score no longer supported for tallies, & &please remove") case ('n1n') - t % score_bins(j) = SCORE_N_1N + if (run_CE) then + t % score_bins(j) = SCORE_N_1N - ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + else + call fatal_error("Cannot tally n1n rate in multi-group mode!") + end if case ('n2n', '(n,2n)') t % score_bins(j) = N_2N @@ -4506,6 +4528,7 @@ contains type(Node), pointer :: doc => null() type(Node), pointer :: node_xsdata => null() type(NodeList), pointer :: node_xsdata_list => null() + real(8), allocatable :: rev_energy_bins(:) ! Check if cross_sections.xml exists inquire(FILE=path_cross_sections, EXIST=file_exists) @@ -4527,6 +4550,7 @@ contains call fatal_error("groups element must exist!") end if + allocate(rev_energy_bins(energy_groups + 1)) allocate(energy_bins(energy_groups + 1)) if (check_for_node(doc, "group_structure")) then ! Get neutron group structure @@ -4535,6 +4559,9 @@ contains call fatal_error("group_structures element must exist!") end if + ! First reverse the order of energy_groups + energy_bins = energy_bins(energy_groups + 1:1:-1) + allocate(energy_bin_avg(energy_groups)) do i = 1, energy_groups energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1)) @@ -4548,7 +4575,7 @@ contains ! If not given, estimate them by using average energy in group which is ! assumed to be the midpoint do i = 1, energy_groups - inverse_velocities(i) = & + inverse_velocities(i) = ONE / & (sqrt(TWO * energy_bin_avg(i) / (MASS_NEUTRON_MEV)) * & C_LIGHT * 100.0_8) end do diff --git a/src/macroxs_header.F90 b/src/macroxs_header.F90 deleted file mode 100644 index 2a12345108..0000000000 --- a/src/macroxs_header.F90 +++ /dev/null @@ -1,860 +0,0 @@ -module macroxs_header - - use constants, only: MAX_FILE_LEN, ZERO, ONE, TWO, PI - use list_header, only: ListInt - use material_header, only: material - use math, only: calc_pn, calc_rn, expand_harmonic, find_angle - use nuclide_header - use random_lcg, only: prn - use scattdata_header - - implicit none - -!=============================================================================== -! MACROXS_* contains cached macroscopic cross sections for the material a -! particle is traveling through -!=============================================================================== - - type, abstract :: MacroXS - ! Data Order - integer :: order - - contains - procedure(macroxs_init_), deferred :: init ! initializes object - procedure(macroxs_get_xs_), deferred :: get_xs ! Return xs - ! Sample the outgoing energy from a fission event - procedure(macroxs_sample_fission_), deferred :: sample_fission_energy - ! Sample the outgoing energy and angle from a scatter event - procedure(macroxs_sample_scatter_), deferred :: sample_scatter - ! Calculate the material specific MGXS data from the nuclides - procedure(macroxs_calculate_xs_), deferred :: calculate_xs - end type MacroXS - - abstract interface - subroutine macroxs_init_(this, mat, nuclides, groups, get_kfiss, get_fiss, & - max_order, scatt_type, legendre_mu_points, & - error_code, error_text) - import MacroXS, Material, NuclideMGContainer, MAX_LINE_LEN - class(MacroXS), intent(inout) :: this ! The MacroXS to initialize - type(Material), pointer, intent(in) :: mat ! base material - type(NuclideMGContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from - integer, intent(in) :: groups ! Number of E groups - logical, intent(in) :: get_kfiss ! Should we get kfiss data? - logical, intent(in) :: get_fiss ! Should we get fiss data? - integer, intent(in) :: max_order ! Maximum requested order - integer, intent(in) :: scatt_type ! Legendre or Tabular Scatt? - integer, intent(in) :: legendre_mu_points ! Treat as Leg or Tabular? - integer, intent(inout) :: error_code ! Code signifying error - character(MAX_LINE_LEN), intent(inout) :: error_text ! Error message to print - end subroutine macroxs_init_ - - function macroxs_get_xs_(this, g, xstype, gout, uvw) result(xs) - import MacroXS - class(MacroXS), intent(in) :: this ! The MacroXS to initialize - integer, intent(in) :: g ! Incoming Energy group - character(*) , intent(in) :: xstype ! Cross Section Type - integer, optional, intent(in) :: gout ! Outgoing Energy group - real(8), optional, intent(in) :: uvw(3) ! Requested Angle - real(8) :: xs ! Resultant xs - end function macroxs_get_xs_ - - function macroxs_sample_fission_(this, gin, uvw) result(gout) - import MacroXS - class(MacroXS), intent(in) :: this ! Data to work with - integer, intent(in) :: gin ! Incoming energy group - real(8), intent(in) :: uvw(3) ! Particle Direction - integer :: gout ! Sampled outgoing group - - end function macroxs_sample_fission_ - - subroutine macroxs_sample_scatter_(this, uvw, gin, gout, mu, wgt) - import MacroXS - class(MacroXS), intent(in) :: this - real(8), intent(in) :: uvw(3) ! Incoming neutron direction - integer, intent(in) :: gin ! Incoming neutron group - integer, intent(out) :: gout ! Sampled outgoin group - real(8), intent(out) :: mu ! Sampled change in angle - real(8), intent(inout) :: wgt ! Particle weight - end subroutine macroxs_sample_scatter_ - - subroutine macroxs_calculate_xs_(this, gin, uvw, xs) - import MacroXS, MaterialMacroXS - class(MacroXS), intent(in) :: this - integer, intent(in) :: gin ! Incoming neutron group - real(8), intent(in) :: uvw(3) ! Incoming neutron direction - type(MaterialMacroXS), intent(inout) :: xs - end subroutine macroxs_calculate_xs_ - end interface - - type, extends(MacroXS) :: MacroXSIso - ! Microscopic cross sections - real(8), allocatable :: total(:) ! total cross section - real(8), allocatable :: absorption(:) ! absorption cross section - class(ScattData), allocatable :: scatter ! scattering information - real(8), allocatable :: nu_fission(:) ! nu-fission - real(8), allocatable :: k_fission(:) ! kappa-fission - real(8), allocatable :: fission(:) ! fission x/s - real(8), allocatable :: scattxs(:) ! scattering xs - real(8), allocatable :: chi(:,:) ! fission spectra - - contains - procedure :: init => macroxsiso_init ! inits object - procedure :: get_xs => macroxsiso_get_xs ! Returns xs - procedure :: sample_fission_energy => macroxsiso_sample_fission_energy - procedure :: sample_scatter => macroxsiso_sample_scatter - procedure :: calculate_xs => macroxsiso_calculate_xs - end type MacroXSIso - - type, extends(MacroXS) :: MacroXSAngle - ! Macroscopic cross sections - real(8), allocatable :: total(:,:,:) ! total cross section - real(8), allocatable :: absorption(:,:,:) ! absorption cross section - type(ScattDataContainer), allocatable :: scatter(:,:) ! scattering information - real(8), allocatable :: nu_fission(:,:,:) ! nu-fission - real(8), allocatable :: k_fission(:,:,:) ! kappa-fission - real(8), allocatable :: fission(:,:,:) ! fission x/s - real(8), allocatable :: chi(:,:,:,:) ! fission spectra - real(8), allocatable :: scattxs(:,:,:) ! scattering xs - real(8), allocatable :: polar(:) ! polar angles - real(8), allocatable :: azimuthal(:) ! azimuthal angles - - contains - procedure :: init => macroxsangle_init ! inits object - procedure :: get_xs => macroxsangle_get_xs ! Returns xs - procedure :: sample_fission_energy => macroxsangle_sample_fission_energy - procedure :: sample_scatter => macroxsangle_sample_scatter - procedure :: calculate_xs => macroxsangle_calculate_xs - end type MacroXSAngle - -!=============================================================================== -! MACROXSCONTAINER pointer array for storing MacroXS objects. -!=============================================================================== - - type MacroXSContainer - class(MacroXS), allocatable :: obj - end type MacroXSContainer - -contains - -!=============================================================================== -! MACROXS*_INIT sets the MacroXS Data -!=============================================================================== - - subroutine macroxsiso_init(this, mat, nuclides, groups, get_kfiss, get_fiss, & - max_order, scatt_type, legendre_mu_points, error_code, error_text) - class(MacroXSIso), intent(inout) :: this ! The MacroXS to initialize - type(Material), pointer, intent(in) :: mat ! base material - type(NuclideMGContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from - integer, intent(in) :: groups ! Number of E groups - logical, intent(in) :: get_kfiss ! Should we get kfiss data? - logical, intent(in) :: get_fiss ! Should we get fiss data? - integer, intent(in) :: max_order ! Maximum requested order - integer, intent(in) :: scatt_type ! How is data presented - integer, intent(in) :: legendre_mu_points ! Treat as Leg or Tabular? - integer, intent(inout) :: error_code ! Code signifying error - character(MAX_LINE_LEN), intent(inout) :: error_text ! Error message to print - - integer :: i ! loop index over nuclides - integer :: gin, gout ! group indices - real(8) :: atom_density ! atom density of a nuclide - integer :: imu - real(8) :: norm - integer :: mat_max_order, order, l - real(8), allocatable :: temp_mult(:,:) - real(8), allocatable :: temp_energy(:,:) - real(8), allocatable :: scatt_coeffs(:,:,:) - - ! Initialize error data - error_code = 0 - error_text = '' - - ! If we have tabular only data, then make sure all datasets have same size - if (scatt_type == ANGLE_HISTOGRAM) then - ! Check all scattering data of same size - order = nuclides(mat % nuclide(1)) % obj % order - do i = 2, mat % n_nuclides - if (order /= nuclides(mat % nuclide(i)) % obj % order) then - error_code = 1 - error_text = "All Histogram Scattering Entries Must Be Same Length!" - return - end if - end do - ! Ok, got our order, store it - this % order = order - - ! Allocate stuff for later - allocate(scatt_coeffs(order, groups, groups)) - scatt_coeffs = ZERO - allocate(ScattDataHistogram :: this % scatter) - - else if (scatt_type == ANGLE_TABULAR) then - ! Check all scattering data of same size - order = nuclides(mat % nuclide(1)) % obj % order - do i = 2, mat % n_nuclides - if (order /= nuclides(mat % nuclide(i)) % obj % order) then - error_code = 1 - error_text = "All Tabular Scattering Entries Must Be Same Length!" - return - end if - end do - ! Ok, got our order, store it - this % order = order - - ! Allocate stuff for later - allocate(scatt_coeffs(order, groups, groups)) - scatt_coeffs = ZERO - allocate(ScattDataTabular :: this % scatter) - - else if (scatt_type == ANGLE_LEGENDRE) then - ! Otherwise find the maximum scattering order - ! Need to determine the maximum scattering order of all data in this material - mat_max_order = 0 - do i = 1, mat % n_nuclides - if (nuclides(mat % nuclide(i)) % obj % order > mat_max_order) then - mat_max_order = nuclides(mat % nuclide(i)) % obj % order - end if - end do - - ! Now need to compare this material maximum scattering order with - ! the problem wide max scatt order and use whichever is lower - order = min(mat_max_order, max_order) - this % order = order + 1 - - ! Now we can allocate our scatt_coeffs object accordingly - allocate(scatt_coeffs(order + 1, groups, groups)) - scatt_coeffs = ZERO - if (legendre_mu_points == 1) then - allocate(ScattDataLegendre :: this % scatter) - else - allocate(ScattDataTabular :: this % scatter) - end if - end if - - ! Allocate and initialize data within macro_xs(i_mat) object - allocate(this % total(groups)) - this % total = ZERO - allocate(this % absorption(groups)) - this % absorption = ZERO - if (get_fiss) then - allocate(this % fission(groups)) - this % fission = ZERO - end if - if (get_kfiss) then - allocate(this % k_fission(groups)) - this % k_fission = ZERO - end if - allocate(this % nu_fission(groups)) - this % nu_fission = ZERO - allocate(this % chi(groups, groups)) - this % chi = ZERO - allocate(temp_energy(groups, groups)) - temp_energy = ZERO - allocate(temp_mult(groups, groups)) - temp_mult = ZERO - allocate(this % scattxs(groups)) - - ! Add contribution from each nuclide in material - do i = 1, mat % n_nuclides - ! Copy atom density of nuclide in material - atom_density = mat % atom_density(i) - - ! Perform our operations which depend upon the type - select type(nuc => nuclides(mat % nuclide(i)) % obj) - type is (NuclideIso) - - ! 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 - if (nuc % fissionable) then - if (allocated(nuc % chi)) then - do gin = 1, groups - do gout = 1, groups - this % chi(gout,gin) = this % chi(gout,gin) + atom_density * & - nuc % chi(gout) * nuc % nu_fission(gin,1) - end do - end do - this % nu_fission = this % nu_fission + atom_density * & - nuc % nu_fission(:,1) - else - this % chi = this % chi + atom_density * nuc % nu_fission - do gin = 1, groups - this % nu_fission(gin) = this % nu_fission(gin) + atom_density * & - sum(nuc % nu_fission(:,gin)) - end do - end if - if (get_fiss) then - this % fission = this % fission + atom_density * nuc % fission - end if - if (get_kfiss) then - this % k_fission = this % k_fission + atom_density * nuc % k_fission - end if - end if - - ! Now time to do the scattering - do gin = 1, groups - do gout = 1, groups - if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then - ! Transfer matrix - temp_energy(gout,gin) = temp_energy(gout,gin) + atom_density * & - sum(nuc % scatter(gout,gin,:)) - - ! Determine the angular distribution - do imu = 1, order - scatt_coeffs(imu, gout, gin) = scatt_coeffs(imu, gout, gin) + & - nuc % scatter(gout,gin,imu) * & - atom_density - end do - - else if (scatt_type == ANGLE_LEGENDRE) then - ! Transfer matrix - temp_energy(gout,gin) = temp_energy(gout,gin) + atom_density * & - nuc % scatter(gout,gin,1) - - ! Determine the angular distribution coefficients so we can later - ! expand do the complete distribution - do l = 1, min(nuc % order, order) + 1 - scatt_coeffs(l, gout, gin) = scatt_coeffs(l, gout, gin) + & - nuc % scatter(gout,gin,l) * & - atom_density - end do - - end if - - ! Multiplicity matrix - temp_mult(gout,gin) = temp_mult(gout,gin) + atom_density * & - nuc % mult(gout,gin) - end do - end do - type is (NuclideAngle) - error_code = 1 - error_text = "Invalid Passing of NuclideAngle to MacroXSIso Object" - return - end select - end do - - ! Store the scattering xs - if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then - this % scattxs(:) = sum(sum(scatt_coeffs(:,:,:),dim=1),dim=1) - else if (scatt_type == ANGLE_LEGENDRE) then - this % scattxs(:) = sum(scatt_coeffs(1,:,:),dim=1) - end if - - ! Normalize the scatt_coeffs - do gin = 1, groups - do gout = 1, groups - if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then - norm = sum(scatt_coeffs(:,gout,gin)) - else if (scatt_type == ANGLE_LEGENDRE) then - norm = scatt_coeffs(1,gout,gin) - end if - if (norm /= ZERO) then - scatt_coeffs(:, gout, gin) = scatt_coeffs(:, gout,gin) / norm - end if - end do - ! Now normalize temp_energy (outgoing scattering energy probabilities) - norm = sum(temp_energy(:,gin)) - if (norm > ZERO) then - temp_energy(:,gin) = temp_energy(:,gin) / norm - end if - end do - - if (scatt_type == ANGLE_LEGENDRE .and. legendre_mu_points /= 1) then - call this % scatter % init(legendre_mu_points, temp_energy, temp_mult, & - scatt_coeffs) - else - call this % scatter % init(this % order, temp_energy, temp_mult, & - scatt_coeffs) - end if - - ! Now normalize chi - if (mat % fissionable) then - do gin = 1, groups - ! Normalize Chi - norm = sum(this % chi(:,gin)) - if (norm > ZERO) then - this % chi(:,gin) = this % chi(:,gin) / norm - end if - end do - end if - - ! Deallocate temporaries for the next material - deallocate(scatt_coeffs, temp_energy, temp_mult) - - end subroutine macroxsiso_init - - subroutine macroxsangle_init(this, mat, nuclides, groups, get_kfiss, get_fiss, & - max_order, scatt_type, legendre_mu_points, error_code, error_text) - class(MacroXSAngle), intent(inout) :: this ! The MacroXS to initialize - type(Material), pointer, intent(in) :: mat ! base material - type(NuclideMGContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from - integer, intent(in) :: groups ! Number of E groups - logical, intent(in) :: get_kfiss ! Should we get kfiss data? - logical, intent(in) :: get_fiss ! Should we get fiss data? - integer, intent(in) :: max_order ! Maximum requested order - integer, intent(in) :: scatt_type ! Legendre or Tabular Scatt? - integer, intent(in) :: legendre_mu_points ! Treat as Leg or Tabular? - integer, intent(inout) :: error_code ! Code signifying error - character(MAX_LINE_LEN), intent(inout) :: error_text ! Error message to print - - integer :: i ! loop index over nuclides - integer :: gin, gout ! group indices - real(8) :: atom_density ! atom density of a nuclide - integer :: ipol, iazi, npol, nazi - integer :: imu - real(8) :: norm - integer :: mat_max_order, order, l - real(8), allocatable :: temp_mult(:,:,:,:) - real(8), allocatable :: temp_energy(:,:,:,:) - real(8), allocatable :: scatt_coeffs(:,:,:,:,:) - - ! Initialize error data - error_code = 0 - error_text = '' - - ! Get the number of each polar and azi angles and make sure all the - ! NuclideAngle types have the same number of these angles - npol = -1 - nazi = -1 - do i = 1, mat % n_nuclides - select type(nuc => nuclides(mat % nuclide(i)) % obj) - type is (NuclideAngle) - if (npol == -1) then - npol = nuc % n_pol - nazi = nuc % n_azi - allocate(this % polar(npol)) - this % polar = nuc % polar - allocate(this % azimuthal(nazi)) - this % azimuthal = nuc % azimuthal - else - if ((npol /= nuc % n_pol) .or. (nazi /= nuc % n_azi)) then - error_code = 1 - error_text = "All Angular Data Must Be Same Length!" - end if - end if - end select - end do - - ! If we have tabular only data, then make sure all datasets have same size - if (scatt_type == ANGLE_HISTOGRAM) then - ! Check all scattering data of same size - order = nuclides(mat % nuclide(1)) % obj % order - do i = 2, mat % n_nuclides - if (order /= nuclides(mat % nuclide(i)) % obj % order) then - error_code = 1 - error_text = "All Histogram Scattering Entries Must Be Same Length!" - return - end if - end do - ! Ok, got our order, store it - this % order = order - - ! Allocate stuff for later - allocate(scatt_coeffs(order, groups, groups, nazi, npol)) - scatt_coeffs = ZERO - allocate(this % scatter(nazi, npol)) - do ipol = 1, npol - do iazi = 1, nazi - allocate(ScattDataHistogram :: this % scatter(iazi, ipol) % obj) - end do - end do - - else if (scatt_type == ANGLE_TABULAR) then - ! Check all scattering data of same size - order = nuclides(mat % nuclide(1)) % obj % order - do i = 2, mat % n_nuclides - if (order /= nuclides(mat % nuclide(i)) % obj % order) then - error_code = 1 - error_text = "All Tabular Scattering Entries Must Be Same Length!" - return - end if - end do - ! Ok, got our order, store it - this % order = order - - ! Allocate stuff for later - allocate(scatt_coeffs(order, groups, groups, nazi, npol)) - scatt_coeffs = ZERO - allocate(this % scatter(nazi, npol)) - do ipol = 1, npol - do iazi = 1, nazi - allocate(ScattDataTabular :: this % scatter(iazi, ipol) % obj) - end do - end do - - else if (scatt_type == ANGLE_LEGENDRE) then - ! Otherwise find the maximum scattering order - ! Need to determine the maximum scattering order of all data in this material - mat_max_order = 0 - do i = 1, mat % n_nuclides - if (nuclides(mat % nuclide(i)) % obj % order > mat_max_order) then - mat_max_order = nuclides(mat % nuclide(i)) % obj % order - end if - end do - - ! Now need to compare this material maximum scattering order with - ! the problem wide max scatt order and use whichever is lower - order = min(mat_max_order, max_order) - this % order = order + 1 - - ! Now we can allocate our scatt_coeffs object accordingly - allocate(scatt_coeffs(order + 1, groups, groups, nazi, npol)) - scatt_coeffs = ZERO - allocate(this % scatter(nazi, npol)) - do ipol = 1, npol - do iazi = 1, nazi - if (legendre_mu_points == 1) then - allocate(ScattDataLegendre :: this % scatter(iazi, ipol) % obj) - else - allocate(ScattDataTabular :: this % scatter(iazi, ipol) % obj) - end if - end do - end do - end if - - ! Allocate and initialize data within macro_xs(i_mat) object - allocate(this % total(groups,nazi,npol)) - this % total = ZERO - allocate(this % absorption(groups,nazi,npol)) - this % absorption = ZERO - if (get_fiss) then - allocate(this % fission(groups,nazi,npol)) - this % fission = ZERO - end if - if (get_kfiss) then - allocate(this % k_fission(groups,nazi,npol)) - this % k_fission = ZERO - end if - allocate(this % nu_fission(groups,nazi,npol)) - this % nu_fission = ZERO - allocate(this % chi(groups, groups, nazi, npol)) - this % chi = ZERO - allocate(temp_energy(groups,groups,nazi,npol)) - temp_energy = ZERO - allocate(temp_mult(groups,groups,nazi,npol)) - temp_mult = ZERO - allocate(this % scattxs(groups,nazi,npol)) - - ! Add contribution from each nuclide in material - do i = 1, mat % n_nuclides - ! Copy atom density of nuclide in material - atom_density = mat % atom_density(i) - - ! Perform our operations which depend upon the type - select type(nuc => nuclides(mat % nuclide(i)) % obj) - type is (NuclideIso) - error_code = 1 - error_text = "Invalid Passing of NuclideIso to MacroXSAngle Object" - return - type is (NuclideAngle) - ! 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 - if (nuc % fissionable) then - if (allocated(nuc % chi)) then - do gin = 1, groups - do gout = 1, groups - this % chi(gout,gin,:,:) = this % chi(gout,gin,:,:) + atom_density * & - nuc % chi(gout,:,:) * nuc % nu_fission(gin,1,:,:) - end do - end do - this % nu_fission = this % nu_fission + atom_density * & - nuc % nu_fission(:,1,:,:) - else - this % chi = this % chi + atom_density * nuc % nu_fission - do gin = 1, groups - this % nu_fission(gin,:,:) = this % nu_fission(gin,:,:) + atom_density * & - sum(nuc % nu_fission(:,gin,:,:),dim=1) - end do - end if - if (get_fiss) then - this % fission = this % fission + atom_density * nuc % fission - end if - if (get_kfiss) then - this % k_fission = this % k_fission + atom_density * nuc % k_fission - end if - end if - - ! Now time to do the scattering - do gin = 1, groups - do gout = 1, groups - if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then - ! Transfer matrix - temp_energy(gout,gin,:,:) = temp_energy(gout,gin,:,:) + atom_density * & - sum(nuc % scatter(gout,gin,:,:,:),dim=1) - - ! Determine the angular distribution - do imu = 1, order - scatt_coeffs(imu,gout,gin,:,:) = scatt_coeffs(imu,gout,gin,:,:) + & - nuc % scatter(gout,gin,imu,:,:) * & - atom_density - end do - else if (scatt_type == ANGLE_LEGENDRE) then - ! Transfer matrix - temp_energy(gout,gin,:,:) = temp_energy(gout,gin,:,:) + atom_density * & - nuc % scatter(gout,gin,1,:,:) - - ! Determine the angular distribution coefficients so we can later - ! expand do the complete distribution - do l = 1, min(nuc % order, order) + 1 - scatt_coeffs(l, gout, gin,:,:) = scatt_coeffs(l, gout, gin,:,:) + & - nuc % scatter(gout,gin,l,:,:) * & - atom_density - end do - end if - - ! Multiplicity matrix - temp_mult(gout,gin,:,:) = temp_mult(gout,gin,:,:) + atom_density * & - nuc % mult(gout,gin,:,:) - end do - end do - end select - end do - - ! Store the scattering xs - if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then - this % scattxs(:,:,:) = sum(sum(scatt_coeffs(:,:,:,:,:),dim=1),dim=1) - else if (scatt_type == ANGLE_LEGENDRE) then - this % scattxs(:,:,:) = sum(scatt_coeffs(1,:,:,:,:),dim=1) - end if - - ! Normalize the scatt_coeffs - do ipol = 1, npol - do iazi = 1, nazi - do gin = 1, groups - do gout = 1, groups - if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then - norm = sum(scatt_coeffs(:,gout,gin,iazi,ipol)) - else if (scatt_type == ANGLE_LEGENDRE) then - norm = scatt_coeffs(1,gout,gin,iazi,ipol) - end if - if (norm /= ZERO) then - scatt_coeffs(:,gout,gin,iazi,ipol) = & - scatt_coeffs(:,gout,gin,iazi,ipol) / norm - end if - end do - ! Now normalize temp_energy (outgoing scattering energy probabilities) - norm = sum(temp_energy(:,gin,iazi,ipol)) - if (norm > ZERO) then - temp_energy(:,gin,iazi,ipol) = temp_energy(:,gin,iazi,ipol) / norm - end if - end do - - if (scatt_type == ANGLE_LEGENDRE .and. legendre_mu_points /= 1) then - call this % scatter(iazi, ipol) % obj % init(legendre_mu_points, & - temp_energy(:,:,iazi,ipol), temp_mult(:,:,iazi,ipol), & - scatt_coeffs(:,:,:,iazi,ipol)) - else - call this % scatter(iazi, ipol) % obj % init(this % order, & - temp_energy(:,:,iazi,ipol), temp_mult(:,:,iazi,ipol), & - scatt_coeffs(:,:,:,iazi,ipol)) - end if - - end do - end do - - ! Now go through and normalize chi - if (mat % fissionable) then - do ipol = 1, npol - do iazi = 1, nazi - do gin = 1, groups - ! Normalize Chi - norm = sum(this % chi(:,gin,iazi,ipol)) - if (norm > ZERO) then - this % chi(:,gin,iazi,ipol) = this % chi(:,gin,iazi,ipol) / norm - end if - end do - end do - end do - end if - - ! Deallocate temporaries for the next material - deallocate(scatt_coeffs, temp_energy, temp_mult) - - end subroutine macroxsangle_init - -!=============================================================================== -! MACROXS_*_GET_XS returns the requested data type -!=============================================================================== - - function macroxsiso_get_xs(this, g, xstype, gout, uvw) result(xs) - class(MacroXSIso), intent(in) :: this ! The MacroXS to initialize - integer, intent(in) :: g ! Incoming Energy group - character(*) , intent(in) :: xstype ! Type of xs requested - integer, optional, intent(in) :: gout ! Outgoing Energy group - real(8), optional, intent(in) :: uvw(3) ! Requested Angle - real(8) :: xs ! Requested x/s - - select case(xstype) - case('total') - xs = this % total(g) - case('absorption') - xs = this % absorption(g) - case('fission') - xs = this % fission(g) - case('k_fission') - xs = this % k_fission(g) - case('nu_fission') - xs = this % nu_fission(g) - case('scatter') - xs = this % scattxs(g) - case('mult') - if (present(gout)) then - xs = this % scatter % mult(gout,g) - else - xs = sum(this % scatter % mult(:,g)) - end if - end select - - end function macroxsiso_get_xs - - function macroxsangle_get_xs(this, g, xstype, gout,uvw) result(xs) - class(MacroXSAngle), intent(in) :: this ! The MacroXS to initialize - integer, intent(in) :: g ! Incoming Energy group - character(*) , intent(in) :: xstype ! Type of xs requested - integer, optional, intent(in) :: gout ! Outgoing Energy group - real(8), optional, intent(in) :: uvw(3) ! Requested Angle - real(8) :: xs ! Requested x/s - - integer :: iazi, ipol - - if (present(uvw)) then - call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) - select case(xstype) - case('total') - xs = this % total(g,iazi,ipol) - case('absorption') - xs = this % absorption(g,iazi,ipol) - case('fission') - xs = this % fission(g,iazi,ipol) - case('k_fission') - xs = this % k_fission(g,iazi,ipol) - case('nu_fission') - xs = this % nu_fission(g,iazi,ipol) - case('scatter') - xs = this % scattxs(g,iazi,ipol) - case('mult') - if (present(gout)) then - xs = this % scatter(iazi,ipol) % obj % mult(gout,g) - else - xs = sum(this % scatter(iazi,ipol) % obj % mult(:,g)) - end if - end select - end if - - end function macroxsangle_get_xs - -!=============================================================================== -! MACROXS_*_SAMPLE_FISSION_ENERGY samples the outgoing energy from a fission -! event -!=============================================================================== - - function macroxsiso_sample_fission_energy(this, gin, uvw) result(gout) - class(MacroXSIso), intent(in) :: this ! Data to work with - integer, intent(in) :: gin ! Incoming energy group - real(8), intent(in) :: uvw(3) ! Particle Direction - integer :: gout ! Sampled outgoing group - real(8) :: xi ! Our random number - real(8) :: prob ! Running probability - - xi = prn() - prob = ZERO - gout = 0 - - do while (prob < xi) - gout = gout + 1 - prob = prob + this % chi(gout,gin) - end do - - end function macroxsiso_sample_fission_energy - - function macroxsangle_sample_fission_energy(this, gin, uvw) result(gout) - class(MacroXSAngle), intent(in) :: this ! Data to work with - integer, intent(in) :: gin ! Incoming energy group - real(8), intent(in) :: uvw(3) ! Particle Direction - integer :: gout ! Sampled outgoing group - real(8) :: xi ! Our random number - real(8) :: prob ! Running probability - integer :: iazi, ipol - - call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) - - xi = prn() - prob = ZERO - gout = 0 - - do while (prob < xi) - gout = gout + 1 - prob = prob + this % chi(gout,gin,iazi,ipol) - end do - - end function macroxsangle_sample_fission_energy - -!=============================================================================== -! MACROXS*_SAMPLE_SCATTER Selects outgoing energy and angle after a scatter -! event -!=============================================================================== - - subroutine macroxsiso_sample_scatter(this, uvw, gin, gout, mu, wgt) - class(MacroXSIso), intent(in) :: this - real(8), intent(in) :: uvw(3) ! Incoming neutron direction - integer, intent(in) :: gin ! Incoming neutron group - integer, intent(out) :: gout ! Sampled outgoin group - real(8), intent(out) :: mu ! Sampled change in angle - real(8), intent(inout) :: wgt ! Particle weight - - call this % scatter % sample(gin, gout, mu, wgt) - - end subroutine macroxsiso_sample_scatter - - subroutine macroxsangle_sample_scatter(this, uvw, gin, gout, mu, wgt) - class(MacroXSAngle), intent(in) :: this - real(8), intent(in) :: uvw(3) ! Incoming neutron direction - integer, intent(in) :: gin ! Incoming neutron group - integer, intent(out) :: gout ! Sampled outgoin group - real(8), intent(out) :: mu ! Sampled change in angle - real(8), intent(inout) :: wgt ! Particle weight - - integer :: iazi, ipol ! Angular indices - - call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) - call this % scatter(iazi,ipol) % obj % sample(gin,gout,mu,wgt) - - end subroutine macroxsangle_sample_scatter - -!=============================================================================== -! MACROXS*_CALCULATE_XS determines the multi-group macroscopic cross sections -! for the material the particle is currently traveling through. -!=============================================================================== - - subroutine macroxsiso_calculate_xs(this, gin, uvw, xs) - class(MacroXSIso), intent(in) :: this - integer, intent(in) :: gin ! Incoming neutron group - real(8), intent(in) :: uvw(3) ! Incoming neutron direction - type(MaterialMacroXS), intent(inout) :: xs ! Resultant MacroXS Data - - xs % total = this % total(gin) - xs % elastic = this % scattxs(gin) - xs % absorption = this % absorption(gin) - xs % nu_fission = this % nu_fission(gin) - - end subroutine macroxsiso_calculate_xs - - subroutine macroxsangle_calculate_xs(this, gin, uvw, xs) - class(MacroXSAngle), intent(in) :: this - integer, intent(in) :: gin ! Incoming neutron group - real(8), intent(in) :: uvw(3) ! Incoming neutron direction - type(MaterialMacroXS), intent(inout) :: xs ! Resultant MacroXS Data - - integer :: iazi, ipol - - call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) - xs % total = this % total(gin, iazi, ipol) - xs % elastic = this % scattxs(gin, iazi, ipol) - xs % absorption = this % absorption(gin, iazi, ipol) - xs % nu_fission = this % nu_fission(gin, iazi, ipol) - - end subroutine macroxsangle_calculate_xs - -end module macroxs_header diff --git a/src/math.F90 b/src/math.F90 index f49220da76..c7ede8d2a5 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -702,31 +702,4 @@ contains end function watt_spectrum - -!=============================================================================== -! find_angle finds the closest angle on the data grid and returns that index -!=============================================================================== - - pure subroutine find_angle(polar, azimuthal, uvw, i_azi, i_pol) - real(8), intent(in) :: polar(:) ! Polar angles [0,pi] - real(8), intent(in) :: azimuthal(:) ! Azi. angles [-pi,pi] - real(8), intent(in) :: uvw(3) ! Direction of motion - integer, intent(inout) :: i_pol ! Closest polar bin - integer, intent(inout) :: i_azi ! Closest azi bin - - real(8) :: my_pol, my_azi, dangle - - ! Convert uvw to polar and azi - - my_pol = acos(uvw(3)) - my_azi = atan2(uvw(2), uvw(1)) - - ! Search for equi-binned angles - dangle = PI / real(size(polar),8) - i_pol = floor(my_pol / dangle + ONE) - dangle = TWO * PI / real(size(azimuthal),8) - i_azi = floor((my_azi + PI) / dangle + ONE) - - end subroutine find_angle - end module math diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 08941870cd..38d6ebc6c5 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -3,9 +3,8 @@ module mgxs_data use constants use error, only: fatal_error use global - use macroxs_header use material_header, only: Material - use nuclide_header + use mgxs_header use output, only: write_message use set_header, only: SetChar use string, only: to_lower @@ -71,7 +70,8 @@ contains if (tallies(i) % score_bins(l) == SCORE_KAPPA_FISSION) then get_kfiss = .true. end if - if (tallies(i) % score_bins(l) == SCORE_FISSION) then + if (tallies(i) % score_bins(l) == SCORE_FISSION .or. & + tallies(i) % score_bins(l) == SCORE_NU_FISSION) then get_fiss = .true. end if end do @@ -118,17 +118,14 @@ contains ! Now allocate accordingly select case(representation) case(MGXS_ISOTROPIC) - allocate(NuclideIso :: nuclides_MG(i_nuclide) % obj) + allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj) case(MGXS_ANGLE) - allocate(NuclideAngle :: nuclides_MG(i_nuclide) % obj) + allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj) end select ! Now read in the data specific to the type we just declared - call nuclides_MG(i_nuclide) % obj % init(node_xsdata, energy_groups, & - get_kfiss, get_fiss) - - ! Keep track of what listing is associated with this nuclide - nuclides_MG(i_nuclide) % obj % listing = i_listing + call nuclides_MG(i_nuclide) % obj % init_file(node_xsdata, & + energy_groups, get_kfiss, get_fiss, max_order, i_listing) ! Add name and alias to dictionary call already_read % add(name) @@ -167,31 +164,8 @@ contains subroutine create_macro_xs() integer :: i_mat ! index in materials array - integer :: i ! loop index over nuclides - integer :: l ! Loop over score bins type(Material), pointer :: mat ! current material - logical :: get_kfiss, get_fiss - integer :: error_code - character(MAX_LINE_LEN) :: error_text integer :: scatt_type - integer :: legendre_mu_points - - ! Find out if we need fission & kappa fission - ! (i.e., are there any SCORE_FISSION or SCORE_KAPPA_FISSION tallies?) - get_kfiss = .false. - get_fiss = .false. - do i = 1, n_tallies - do l = 1, tallies(i) % n_score_bins - if (tallies(i) % score_bins(l) == SCORE_KAPPA_FISSION) then - get_kfiss = .true. - end if - if (tallies(i) % score_bins(l) == SCORE_FISSION) then - get_fiss = .true. - end if - end do - if (get_kfiss .and. get_fiss) & - exit - end do allocate(macro_xs(n_materials)) @@ -203,21 +177,15 @@ contains ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs) ! At the same time, we will find the scattering type, as that will dictate ! how we allocate the scatter object within macroxs - legendre_mu_points = nuclides_MG(mat % nuclide(1)) % obj % legendre_mu_points scatt_type = nuclides_MG(mat % nuclide(1)) % obj % scatt_type select type(nuc => nuclides_MG(mat % nuclide(1)) % obj) - type is (NuclideIso) - allocate(MacroXSIso :: macro_xs(i_mat) % obj) - type is (NuclideAngle) - allocate(MacroXSAngle :: macro_xs(i_mat) % obj) + type is (MgxsIso) + allocate(MgxsIso :: macro_xs(i_mat) % obj) + type is (MgxsAngle) + allocate(MgxsAngle :: macro_xs(i_mat) % obj) end select - - call macro_xs(i_mat) % obj % init(mat, nuclides_MG, energy_groups, & - get_kfiss, get_fiss, max_order, & - scatt_type, legendre_mu_points, & - error_code, error_text) - ! Handle any errors - if (error_code /= 0) call fatal_error(trim(error_text)) + call macro_xs(i_mat) % obj % combine(mat, nuclides_MG, energy_groups, & + max_order, scatt_type, i_mat) end do end subroutine create_macro_xs diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 new file mode 100644 index 0000000000..88c1b23e25 --- /dev/null +++ b/src/mgxs_header.F90 @@ -0,0 +1,1929 @@ +module mgxs_header + + use constants, only: MAX_FILE_LEN, ZERO, ONE, TWO, PI + use error, only: fatal_error + use, intrinsic :: ISO_FORTRAN_ENV, only: OUTPUT_UNIT + use list_header, only: ListInt + use material_header, only: material + use math, only: calc_pn, calc_rn, expand_harmonic, & + evaluate_legendre + use nuclide_header, only: MaterialMacroXS + use random_lcg, only: prn + use scattdata_header + use string + use xml_interface + +!=============================================================================== +! MGXS contains the base mgxs data for a nuclide/material +!=============================================================================== + + type, abstract :: Mgxs + character(len=104) :: name ! name of dataset, e.g. 92235.03c + integer :: zaid ! Z and A identifier, e.g. 92235 + real(8) :: awr ! Atomic Weight Ratio + integer :: listing ! index in xs_listings + real(8) :: kT ! temperature in MeV (k*T) + + ! Fission information + logical :: fissionable ! mgxs object is fissionable? + integer :: scatt_type ! either legendre, histogram, or tabular. + + contains + procedure(mgxs_init_file_), deferred :: init_file ! Initialize the data + procedure(mgxs_print_), deferred :: print ! Writes object info + procedure(mgxs_get_xs_), deferred :: get_xs ! Get the requested xs + procedure(mgxs_combine_), deferred :: combine ! initializes object + ! Sample the outgoing energy from a fission event + procedure(mgxs_sample_fission_), deferred :: sample_fission_energy + ! Sample the outgoing energy and angle from a scatter event + procedure(mgxs_sample_scatter_), deferred :: sample_scatter + ! Calculate the material specific MGXS data from the nuclides + procedure(mgxs_calculate_xs_), deferred :: calculate_xs + end type Mgxs + +!=============================================================================== +! MGXSCONTAINER pointer array for storing Nuclides +!=============================================================================== + + type MgxsContainer + class(Mgxs), pointer :: obj + end type MgxsContainer + +!=============================================================================== +! Interfaces for MGXS +!=============================================================================== + + abstract interface + subroutine mgxs_init_file_(this,node_xsdata,groups,get_kfiss,get_fiss, & + max_order,i_listing) + import Mgxs, Node + class(Mgxs), intent(inout) :: this ! Working Object + type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml + integer, intent(in) :: groups ! Number of Energy groups + logical, intent(in) :: get_kfiss ! Need Kappa-Fission? + logical, intent(in) :: get_fiss ! Should we get fiss data? + integer, intent(in) :: max_order ! Maximum requested order + integer, intent(in) :: i_listing ! Index of listings array + end subroutine mgxs_init_file_ + + subroutine mgxs_print_(this, unit) + import Mgxs + class(Mgxs),intent(in) :: this + integer, optional, intent(in) :: unit + end subroutine mgxs_print_ + + pure function mgxs_get_xs_(this,xstype,gin,gout,uvw,mu) result(xs) + import Mgxs + class(Mgxs), intent(in) :: this + character(*), intent(in) :: xstype ! Cross Section Type + integer, intent(in) :: gin ! Incoming Energy group + integer, optional, intent(in) :: gout ! Outgoing Group + real(8), optional, intent(in) :: uvw(3) ! Requested Angle + real(8), optional, intent(in) :: mu ! Change in angle + real(8) :: xs ! Resultant xs + end function mgxs_get_xs_ + + pure function mgxs_calc_f_(this,gin,gout,mu,uvw,iazi,ipol) result(f) + import Mgxs + class(Mgxs), intent(in) :: this + integer, intent(in) :: gin ! Incoming Energy Group + integer, intent(in) :: gout ! Outgoing Energy Group + real(8), intent(in) :: mu ! Angle of interest + real(8), intent(in), optional :: uvw(3) ! Direction vector + integer, intent(in), optional :: iazi ! Incoming Energy Group + integer, intent(in), optional :: ipol ! Outgoing Energy Group + real(8) :: f ! Return value of f(mu) + + end function mgxs_calc_f_ + + subroutine mgxs_combine_(this,mat,nuclides,groups,max_order,scatt_type, & + i_listing) + import Mgxs, Material, MgxsContainer + class(Mgxs), intent(inout) :: this ! The Mgxs to initialize + type(Material), pointer, intent(in) :: mat ! base material + type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from + integer, intent(in) :: groups ! Number of E groups + integer, intent(in) :: max_order ! Maximum requested order + integer, intent(in) :: scatt_type ! Legendre or Tabular Scatt? + integer, intent(in) :: i_listing ! Index in listings + end subroutine mgxs_combine_ + + function mgxs_sample_fission_(this, gin, uvw) result(gout) + import Mgxs + class(Mgxs), intent(in) :: this ! Data to work with + integer, intent(in) :: gin ! Incoming energy group + real(8), intent(in) :: uvw(3) ! Particle Direction + integer :: gout ! Sampled outgoing group + + end function mgxs_sample_fission_ + + subroutine mgxs_sample_scatter_(this, uvw, gin, gout, mu, wgt) + import Mgxs + class(Mgxs), intent(in) :: this + real(8), intent(in) :: uvw(3) ! Incoming neutron direction + integer, intent(in) :: gin ! Incoming neutron group + integer, intent(out) :: gout ! Sampled outgoin group + real(8), intent(out) :: mu ! Sampled change in angle + real(8), intent(inout) :: wgt ! Particle weight + end subroutine mgxs_sample_scatter_ + + subroutine mgxs_calculate_xs_(this, gin, uvw, xs) + import Mgxs, MaterialMacroXS + class(Mgxs), intent(in) :: this + integer, intent(in) :: gin ! Incoming neutron group + real(8), intent(in) :: uvw(3) ! Incoming neutron direction + type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data + end subroutine mgxs_calculate_xs_ + end interface + +!=============================================================================== +! MGXSISO contains the base MGXS data specifically for +! isotropically weighted MGXS +!=============================================================================== + + type, extends(Mgxs) :: MgxsIso + + ! Microscopic cross sections + real(8), allocatable :: total(:) ! total cross section + real(8), allocatable :: absorption(:) ! absorption cross section + class(ScattData), allocatable :: scatter ! scattering information + real(8), allocatable :: nu_fission(:) ! fission matrix (Gout x Gin) + real(8), allocatable :: k_fission(:) ! kappa-fission + real(8), allocatable :: fission(:) ! neutron production + real(8), allocatable :: chi(:, :) ! Fission Spectra + + contains + procedure :: init_file => mgxsiso_init_file ! Initialize Nuclidic MGXS Data + procedure :: print => mgxsiso_print ! Writes nuclide info + procedure :: get_xs => mgxsiso_get_xs ! Gets Size of Data w/in Object + procedure :: combine => mgxsiso_combine ! inits object + procedure :: sample_fission_energy => mgxsiso_sample_fission_energy + procedure :: sample_scatter => mgxsiso_sample_scatter + procedure :: calculate_xs => mgxsiso_calculate_xs + end type MgxsIso + +!=============================================================================== +! MGXSANGLE contains the base MGXS data specifically for +! angular flux weighted MGXS +!=============================================================================== + + type, extends(Mgxs) :: MgxsAngle + + ! Microscopic cross sections + real(8), allocatable :: total(:, :, :) ! total cross section + real(8), allocatable :: absorption(:, :, :) ! absorption cross section + type(ScattDataContainer), allocatable :: scatter(:, :) ! scattering information + real(8), allocatable :: nu_fission(:, :, :) ! fission matrix (Gout x Gin) + real(8), allocatable :: k_fission(:, :, :) ! kappa-fission + real(8), allocatable :: fission(:, :, :) ! neutron production + real(8), allocatable :: chi(:, :, :, :) ! Fission Spectra + ! In all cases, right-most indices are theta, phi + integer :: n_pol ! Number of polar angles + integer :: n_azi ! Number of azimuthal angles + real(8), allocatable :: polar(:) ! polar angles + real(8), allocatable :: azimuthal(:) ! azimuthal angles + + contains + procedure :: init_file => mgxsang_init_file ! Initialize Nuclidic MGXS Data + procedure :: print => mgxsang_print ! Writes nuclide info + procedure :: get_xs => mgxsang_get_xs ! Gets Size of Data w/in Object + procedure :: combine => mgxsang_combine ! inits object + procedure :: sample_fission_energy => mgxsang_sample_fission_energy + procedure :: sample_scatter => mgxsang_sample_scatter + procedure :: calculate_xs => mgxsang_calculate_xs + end type MgxsAngle + + contains + +!=============================================================================== +! MGXS*_INIT reads in the data from the XML file. At the point of entry +! the file would have been opened and metadata read. This routine begins with +! the xsdata object node itself. +!=============================================================================== + + 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 + + character(MAX_LINE_LEN) :: temp_str + + ! Load the nuclide metadata + call get_node_value(node_xsdata, "name", this % name) + this % name = to_lower(this % name) + if (check_for_node(node_xsdata, "kT")) then + call get_node_value(node_xsdata, "kT", this % kT) + else + this % kT = ZERO + end if + if (check_for_node(node_xsdata, "zaid")) then + call get_node_value(node_xsdata, "zaid", this % zaid) + else + this % zaid = 0 + end if + if (check_for_node(node_xsdata, "awr")) then + call get_node_value(node_xsdata, "awr", this % awr) + else + this % awr = -ONE + end if + if (check_for_node(node_xsdata, "scatt_type")) then + call get_node_value(node_xsdata, "scatt_type", temp_str) + temp_str = trim(to_lower(temp_str)) + if (temp_str == 'legendre') then + this % scatt_type = ANGLE_LEGENDRE + else if (temp_str == 'histogram') then + this % scatt_type = ANGLE_HISTOGRAM + else if (temp_str == 'tabular') then + this % scatt_type = ANGLE_TABULAR + else + call fatal_error("Invalid scatt_type option!") + end if + else + this % scatt_type = ANGLE_LEGENDRE + end if + + if (check_for_node(node_xsdata, "fissionable")) then + call get_node_value(node_xsdata, "fissionable", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then + this % fissionable = .true. + else + this % fissionable = .false. + end if + else + call fatal_error("Fissionable element must be set!") + end if + + ! Keep track of what listing is associated with this nuclide + this % listing = i_listing + + end subroutine mgxs_init_file + + 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 + logical, intent(in) :: get_kfiss ! Need Kappa-Fission? + logical, intent(in) :: get_fiss ! Need fiss data? + integer, intent(in) :: max_order ! Maximum requested order + integer, intent(in) :: i_listing ! Index in listings array + + type(Node), pointer :: node_legendre_mu + character(MAX_LINE_LEN) :: temp_str + logical :: enable_leg_mu + 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, arr_len + integer :: legendre_mu_points, imu + + ! Call generic data gathering routine (will populate the metadata) + call mgxs_init_file(this, node_xsdata, i_listing) + + ! Load the more specific data + allocate(this % nu_fission(groups)) + allocate(this % chi(groups,groups)) + if (this % fissionable) then + if (check_for_node(node_xsdata, "chi")) then + ! Chi was provided, that means they are giving chi and nu-fission + ! vectors + ! Get chi + allocate(temp_arr(groups)) + call get_node_array(node_xsdata, "chi", temp_arr) + do gin = 1, groups + do gout = 1, groups + this % chi(gout, gin) = temp_arr(gout) + end do + ! Normalize chi so its CDF goes to 1 + this % chi(:, gin) = this % chi(:, gin) / sum(this % chi(:, gin)) + end do + deallocate(temp_arr) + + ! Get nu_fission (as a vector) + if (check_for_node(node_xsdata, "nu_fission")) then + call get_node_array(node_xsdata, "nu_fission", this % nu_fission) + else + call fatal_error("If fissionable, must provide nu_fission!") + end if + + else + ! chi isnt provided but is within nu_fission, existing as a matrix + ! So, get nu_fission (as a matrix) + if (check_for_node(node_xsdata, "nu_fission")) then + allocate(temp_arr(groups*groups)) + call get_node_array(node_xsdata, "nu_fission", temp_arr) + allocate(temp_2d(groups, groups)) + temp_2d = reshape(temp_arr, (/groups, groups/)) + deallocate(temp_arr) + else + call fatal_error("If fissionable, must provide nu_fission!") + end if + + ! Set the vector nu-fission from the matrix nu-fission + do gin = 1, groups + this % nu_fission(gin) = sum(temp_2d(:, gin)) + end do + + ! Now pull out information needed for chi + 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)) + end do + deallocate(temp_2d) + end if + ! If we have a need* for the fission and kappa-fission x/s, get them + ! (*Need is defined as will be using it to tally) + if (get_fiss) then + allocate(this % fission(groups)) + if (check_for_node(node_xsdata, "fission")) then + call get_node_array(node_xsdata, "fission", this % fission) + else + call fatal_error("Fission data missing, required due to fission& + & tallies in tallies.xml file!") + end if + end if + if (get_kfiss) then + allocate(this % k_fission(groups)) + if (check_for_node(node_xsdata, "kappa_fission")) then + call get_node_array(node_xsdata, "kappa_fission", this % k_fission) + else + call fatal_error("kappa_fission data missing, required due to & + &kappa-fission tallies in tallies.xml file!") + end if + end if + else + this % nu_fission = ZERO + this % chi = ZERO + end if + + allocate(this % absorption(groups)) + if (check_for_node(node_xsdata, "absorption")) then + call get_node_array(node_xsdata, "absorption", this % absorption) + else + call fatal_error("Must provide absorption!") + end if + + ! Get multiplication data if present + allocate(temp_mult(groups, groups)) + if (check_for_node(node_xsdata, "multiplicity")) then + arr_len = get_arraysize_double(node_xsdata, "multiplicity") + 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/)) + deallocate(temp_arr) + else + call fatal_error("Multiplicity length not same as number of groups& + & squared!") + end if + else + temp_mult(:, :) = ONE + end if + + ! Get scattering treatment information + ! Tabular_legendre tells us if we are to treat the provided + ! Legendre polynomials as tabular data (if enable is true) or leaving + ! them as Legendres (if enable is false, or the default) + + ! Set the default (leave as Legendre polynomials) + enable_leg_mu = .false. + if (check_for_node(node_xsdata, "tabular_legendre")) then + call get_node_ptr(node_xsdata, "tabular_legendre", node_legendre_mu) + if (check_for_node(node_legendre_mu, "enable")) then + call get_node_value(node_legendre_mu, "enable", temp_str) + temp_str = trim(to_lower(temp_str)) + if (temp_str == 'true' .or. temp_str == '1') then + enable_leg_mu = .true. + elseif (temp_str == 'false' .or. temp_str == '0') then + enable_leg_mu = .false. + else + call fatal_error("Unrecognized tabular_legendre/enable: " & + // temp_str) + end if + end if + ! Ok, so if we need to convert to a tabular form, get the user provided + ! number of points + if (enable_leg_mu) then + if (check_for_node(node_legendre_mu, "num_points")) then + call get_node_value(node_legendre_mu, "num_points", & + legendre_mu_points) + if (legendre_mu_points <= 0) & + call fatal_error("num_points element must be positive& + & and non-zero!") + else + ! Set the default number of points (0.0625 spacing) + legendre_mu_points = 33 + end if + end if + end if + + ! Get the library's value for the order + if (check_for_node(node_xsdata, "order")) then + call get_node_value(node_xsdata, "order", order) + else + call fatal_error("Order must be provided!") + end if + + ! Before retrieving the data, store the dimensionality of the data in + ! order_dim. For Legendre data, we usually refer to it as Pn where + ! n is the order. However Pn has n+1 sets of points (since you need to + ! the count the P0 moment). Adjust for that. Histogram and Tabular + ! formats dont need this adjustment. + if (this % scatt_type == ANGLE_LEGENDRE) then + order_dim = order + 1 + else + order_dim = order + end if + + ! 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_for_node(node_xsdata, "scatter")) then + allocate(temp_arr(groups * groups * order_dim)) + call get_node_array(node_xsdata, "scatter", temp_arr) + input_scatt = reshape(temp_arr, (/groups, groups, order_dim/)) + 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 % scatt_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 this is also + ! the default behavior), convert that now. + if (this % scatt_type == ANGLE_LEGENDRE .and. enable_leg_mu) then + ! Convert input parameters to what we need for the rest. + this % scatt_type = ANGLE_TABULAR + order_dim = legendre_mu_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 + ! 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 + end do + 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) + else + call fatal_error("Must provide scatter!") + end if + + ! Allocate and initialize our ScattData Object. + if (this % scatt_type == ANGLE_HISTOGRAM) then + allocate(ScattDataHistogram :: this % scatter) + else if (this % scatt_type == ANGLE_TABULAR) then + allocate(ScattDataTabular :: this % scatter) + else if (this % scatt_type == ANGLE_LEGENDRE) then + allocate(ScattDataLegendre :: this % scatter) + end if + + ! Initialize the ScattData Object + call this % scatter % init(temp_mult, scatt_coeffs) + + ! 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 gin = 1, groups + if (this % absorption(gin) == ZERO) this % absorption(gin) = 1E-10_8 + end do + + ! Get, or infer, total xs data. + allocate(this % total(groups)) + 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(:) + end if + + ! Deallocate temporaries for the next material + 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 + do gin = 1, groups + if (this % total(gin) == ZERO) this % total(gin) = 1E-10_8 + end do + + + end subroutine mgxsiso_init_file + + 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 + logical, intent(in) :: get_kfiss ! Need Kappa-Fission? + logical, intent(in) :: get_fiss ! Should we get fiss data? + integer, intent(in) :: max_order ! Maximum requested order + integer, intent(in) :: i_listing ! Index in listings array + + type(Node), pointer :: node_legendre_mu + character(MAX_LINE_LEN) :: temp_str + logical :: enable_leg_mu + 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, dangle + integer :: order, order_dim, gin, gout, l, arr_len + 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) + + 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!") + 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!") + end if + + ! Load angle data, if present (else equally spaced) + allocate(this % polar(this % n_pol)) + allocate(this % azimuthal(this % n_azi)) + if (check_for_node(node_xsdata, "polar")) then + call fatal_error("User-Specified polar angle bins not yet supported!") + ! When this feature is supported, this line will be activated + call get_node_array(node_xsdata, "polar", this % polar) + else + dangle = PI / real(this % n_pol, 8) + do ipol = 1, this % n_pol + this % polar(ipol) = (real(ipol, 8) - HALF) * dangle + end do + end if + if (check_for_node(node_xsdata, "azimuthal")) then + call fatal_error("User-Specified azimuthal angle bins not yet supported!") + ! When this feature is supported, this line will be activated + call get_node_array(node_xsdata, "azimuthal", this % azimuthal) + else + dangle = TWO * PI / real(this % n_azi, 8) + do iazi = 1, this % n_azi + this % azimuthal(iazi) = -PI + (real(iazi, 8) - HALF) * dangle + end do + end if + + ! Load the more specific data + allocate(this % nu_fission(groups, this % n_azi, this % n_pol)) + allocate(this % chi(groups, groups, this % n_azi, this % n_pol)) + if (this % fissionable) then + if (check_for_node(node_xsdata, "chi")) then + ! Chi was provided, that means they are giving chi and nu-fission + ! vectors + ! Get chi + allocate(temp_arr(1 * groups * this % n_azi * this % n_pol)) + call get_node_array(node_xsdata, "chi", temp_arr) + ! Initialize counter for temp_arr + l = 0 + gin = 1 + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gout = 1, groups + l = l + 1 + this % chi(gout, gin, iazi, ipol) = temp_arr(l) + end do + ! Normalize chi so its CDF goes to 1 + this % chi(:, gin, iazi, ipol) = & + this % chi(:, gin, iazi, ipol) / & + sum(this % chi(:, gin, iazi, ipol)) + end do + end do + + ! Now set all the other gin values + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 2, groups + this % chi(:, gin, iazi, ipol) = & + this % chi(:, 1, iazi, ipol) + end do + end do + end do + deallocate(temp_arr) + + ! Get nu_fission (as a vector) + 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/)) + deallocate(temp_arr) + else + call fatal_error("If fissionable, must provide nu_fission!") + end if + + else + ! chi isnt provided but is within nu_fission, existing as a matrix + ! So, get nu_fission (as a matrix) + if (check_for_node(node_xsdata, "nu_fission")) then + 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/)) + deallocate(temp_arr) + else + call fatal_error("If fissionable, must provide nu_fission!") + end if + + ! Set the vector nu-fission from the matrix nu-fission + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, groups + this % nu_fission(gin, iazi, ipol) = & + sum(temp_4d(:, gin, iazi, ipol)) + end do + end do + end do + + ! Now pull out information needed for chi + this % chi = temp_4d + ! Normalize chi so its CDF goes to 1 + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, groups + this % chi(:, gin, iazi, ipol) = & + this % chi(:, gin, iazi, ipol) / & + sum(this % chi(:, gin, iazi, ipol)) + end do + end do + end do + deallocate(temp_4d) + end if + + ! If we have a need* for the fission and kappa-fission x/s, get them + ! (*Need is defined as will be using it to tally) + if (get_fiss) then + if (check_for_node(node_xsdata, "fission")) then + 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/)) + deallocate(temp_arr) + else + call fatal_error("Fission data missing, required due to fission& + & tallies in tallies.xml file!") + end if + end if + if (get_kfiss) then + if (check_for_node(node_xsdata, "kappa_fission")) then + 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/)) + deallocate(temp_arr) + else + call fatal_error("kappa_fission data missing, required due to & + &kappa-fission tallies in tallies.xml file!") + end if + end if + else + 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/)) + deallocate(temp_arr) + else + 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_for_node(node_xsdata, "multiplicity")) then + arr_len = get_arraysize_double(node_xsdata, "multiplicity") + 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/)) + deallocate(temp_arr) + else + call fatal_error("Multiplicity length not same as number of groups& + & squared!") + end if + else + temp_mult(:, :, :, :) = ONE + end if + + ! Get scattering treatment information + ! Tabular_legendre tells us if we are to treat the provided + ! Legendre polynomials as tabular data (if enable is true) or leaving + ! them as Legendres (if enable is false, or the default) + + ! Set the default (leave as Legendre polynomials) + enable_leg_mu = .false. + if (check_for_node(node_xsdata, "tabular_legendre")) then + call get_node_ptr(node_xsdata, "tabular_legendre", node_legendre_mu) + if (check_for_node(node_legendre_mu, "enable")) then + call get_node_value(node_legendre_mu, "enable", temp_str) + temp_str = trim(to_lower(temp_str)) + if (temp_str == 'true' .or. temp_str == '1') then + enable_leg_mu = .true. + elseif (temp_str == 'false' .or. temp_str == '0') then + enable_leg_mu = .false. + else + call fatal_error("Unrecognized tabular_legendre/enable: " & + // temp_str) + end if + end if + ! Ok, so if we need to convert to a tabular form, get the user provided + ! number of points + if (enable_leg_mu) then + if (check_for_node(node_legendre_mu, "num_points")) then + call get_node_value(node_legendre_mu, "num_points", & + legendre_mu_points) + if (legendre_mu_points <= 0) & + call fatal_error("num_points element must be positive& + & and non-zero!") + else + ! Set the default number of points (0.0625 spacing) + legendre_mu_points = 33 + end if + end if + end if + + ! Get the library's value for the order + if (check_for_node(node_xsdata, "order")) then + call get_node_value(node_xsdata, "order", order) + else + call fatal_error("Order must be provided!") + end if + + ! Before retrieving the data, store the dimensionality of the data in + ! order_dim. For Legendre data, we usually refer to it as Pn where + ! n is the order. However Pn has n+1 sets of points (since you need to + ! the count the P0 moment). Adjust for that. Histogram and Tabular + ! formats dont need this adjustment. + if (this % scatt_type == ANGLE_LEGENDRE) then + order_dim = order + 1 + else + order_dim = order + end if + + ! 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_for_node(node_xsdata, "scatter")) then + 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/)) + 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 % scatt_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 this is also + ! the default behavior), convert that now. + if (this % scatt_type == ANGLE_LEGENDRE .and. enable_leg_mu) then + + ! Convert input parameters to what we need for the rest. + this % scatt_type = ANGLE_TABULAR + order_dim = legendre_mu_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) + else + call fatal_error("Must provide scatter!") + end if + + allocate(this % scatter(this % n_azi, this % n_pol)) + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + ! Allocate and initialize our ScattData Object. + if (this % scatt_type == ANGLE_HISTOGRAM) then + allocate(ScattDataHistogram :: this % scatter(iazi, ipol) % obj) + else if (this % scatt_type == ANGLE_TABULAR) then + allocate(ScattDataTabular :: this % scatter(iazi, ipol) % obj) + else if (this % scatt_type == ANGLE_LEGENDRE) then + allocate(ScattDataLegendre :: this % scatter(iazi, ipol) % obj) + end if + + ! Initialize the ScattData Object + call this % scatter(iazi, ipol) % obj % init(& + temp_mult(:, :, iazi, ipol), & + scatt_coeffs(:, :, :, iazi, ipol)) + end do + end do + ! Deallocate temporaries for the next material + deallocate(input_scatt, scatt_coeffs, temp_mult) + + allocate(this % total(groups, this % n_azi, this % n_pol)) + 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/)) + deallocate(temp_arr) + else + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + this % total(:, iazi, ipol) = this % absorption(:, iazi, ipol) + & + this % scatter(iazi, ipol) % obj % scattxs(:) + end do + end do + end if + + end subroutine mgxsang_init_file + +!=============================================================================== +! MGXS*_PRINT displays information about a continuous-energy neutron +! cross_section table and its reactions and secondary angle/energy distributions +!=============================================================================== + + subroutine mgxs_print(this, unit_) + class(Mgxs), intent(in) :: this + integer, intent(in) :: unit_ + + character(MAX_LINE_LEN) :: temp_str + + ! Basic nuclide information + write(unit_,*) 'MGXS Entry: ' // trim(this % name) + if (this % zaid > 0) then + write(unit_,*) ' ZAID = ' // trim(to_str(this % zaid)) + else if (this % zaid < 0) then + write(unit_,*) ' Material id = ' // trim(to_str(-this % zaid)) + end if + if (this % awr > ZERO) then + write(unit_,*) ' AWR = ' // trim(to_str(this % awr)) + end if + if (this % kT > ZERO) then + write(unit_,*) ' kT = ' // trim(to_str(this % kT)) + end if + if (this % scatt_type == ANGLE_LEGENDRE) then + temp_str = "Legendre" + write(unit_,*) ' Scattering Type = ' // trim(temp_str) + select type(this) + type is (MgxsIso) + temp_str = to_str(size(this % scatter % dist(1) % data,dim=1) - 1) + end select + write(unit_,*) ' Scattering Order = ' // trim(temp_str) + else if (this % scatt_type == ANGLE_HISTOGRAM) then + temp_str = "Histogram" + write(unit_,*) ' Scattering Type = ' // trim(temp_str) + select type(this) + type is (MgxsIso) + temp_str = to_str(size(this % scatter % dist(1) % data,dim=1)) + end select + write(unit_,*) ' Num. Distribution Bins = ' // trim(temp_str) + else if (this % scatt_type == ANGLE_TABULAR) then + temp_str = "Tabular" + write(unit_,*) ' Scattering Type = ' // trim(temp_str) + select type(this) + type is (MgxsIso) + temp_str = to_str(size(this % scatter % dist(1) % data,dim=1)) + end select + write(unit_,*) ' Num. Distribution Points = ' // trim(temp_str) + end if + write(unit_,*) ' Fissionable = ', this % fissionable + + end subroutine mgxs_print + + subroutine mgxsiso_print(this, unit) + + class(MgxsIso), intent(in) :: this + integer, optional, intent(in) :: unit + + integer :: unit_ ! unit to write to + integer :: size_total, size_scattmat, size_mgxs + integer :: gin + + ! set default unit for writing information + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + ! Write Basic Nuclide Information + call mgxs_print(this, unit_) + + ! Determine size of mgxs and scattering matrices + size_scattmat = 0 + do gin = 1, size(this % scatter % energy) + size_scattmat = size_scattmat + & + 2 * size(this % scatter % energy(gin) % data) + & + size(this % scatter % dist(gin) % data) + end do + size_scattmat = size_scattmat + size(this % scatter % scattxs) + size_scattmat = size_scattmat * 8 + + size_mgxs = size(this % total) + size(this % absorption) + & + size(this % nu_fission) + size(this % k_fission) + & + size(this % fission) + size(this % chi) + size_mgxs = size_mgxs * 8 + + ! Calculate total memory + size_total = size_scattmat + size_mgxs + + ! Write memory used + write(unit_,*) ' Memory Requirements' + write(unit_,*) ' Cross sections = ' // trim(to_str(size_mgxs)) // ' bytes' + write(unit_,*) ' Scattering Matrices = ' // & + trim(to_str(size_scattmat)) // ' bytes' + write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' + + ! Blank line at end of nuclide + write(unit_,*) + + end subroutine mgxsiso_print + + subroutine mgxsang_print(this, unit) + + class(MgxsAngle), intent(in) :: this + integer, optional, intent(in) :: unit + + integer :: unit_ ! unit to write to + integer :: size_total, size_scattmat, size_mgxs + integer :: ipol, iazi, gin + + ! set default unit for writing information + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + ! Write Basic Nuclide Information + call mgxs_print(this, unit_) + + write(unit_,*) ' # of Polar Angles = ' // trim(to_str(this % n_pol)) + write(unit_,*) ' # of Azimuthal Angles = ' // trim(to_str(this % n_azi)) + + ! Determine size of mgxs and scattering matrices + size_scattmat = 0 + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, size(this % scatter(iazi, ipol) % obj % energy) + size_scattmat = size_scattmat + & + 2 * size(this % scatter(iazi, ipol) % obj % energy(gin) % data) + & + size(this % scatter(iazi, ipol) % obj % dist(gin) % data) + end do + size_scattmat = size_scattmat + & + size(this % scatter(iazi, ipol) % obj % scattxs) + end do + end do + size_scattmat = size_scattmat * 8 + + size_mgxs = size(this % total) + size(this % absorption) + & + size(this % nu_fission) + size(this % k_fission) + & + size(this % fission) + size(this % chi) + size_mgxs = size_mgxs * 8 + + ! Calculate total memory + size_total = size_scattmat + size_mgxs + + ! Write memory used + write(unit_,*) ' Memory Requirements' + write(unit_,*) ' Cross sections = ' // trim(to_str(size_mgxs)) // ' bytes' + write(unit_,*) ' Scattering Matrices = ' // & + trim(to_str(size_scattmat)) // ' bytes' + write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' + + ! Blank line at end of nuclide + write(unit_,*) + end subroutine mgxsang_print + +!=============================================================================== +! MGXS*_GET_XS returns the requested data cross section data +!=============================================================================== + + pure function mgxsiso_get_xs(this, xstype, gin, gout, uvw, mu) result(xs) + class(MgxsIso), intent(in) :: this ! The Mgxs to initialize + character(*) , intent(in) :: xstype ! Type of xs requested + integer, intent(in) :: gin ! Incoming Energy group + integer, optional, intent(in) :: gout ! Outgoing Energy group + real(8), optional, intent(in) :: uvw(3) ! Requested Angle + real(8), optional, intent(in) :: mu ! Change in angle + real(8) :: xs ! Requested x/s + + select case(xstype) + case('total') + xs = this % total(gin) + case('absorption') + xs = this % absorption(gin) + case('fission') + if (allocated(this % fission)) then + xs = this % fission(gin) + else + xs = ZERO + end if + case('kappa_fission') + if (allocated(this % k_fission)) then + xs = this % k_fission(gin) + else + xs = ZERO + end if + case('nu_fission') + xs = this % nu_fission(gin) + case('chi') + if (present(gout)) then + xs = this % chi(gout,gin) + else + ! Not sure youd want a 1 or a 0, but here you go! + xs = sum(this % chi(:, gin)) + end if + case('scatter') + if (present(gout)) then + if (gout < this % scatter % gmin(gin) .or. & + gout > this % scatter % gmax(gin)) then + xs = ZERO + else + xs = this % scatter % scattxs(gin) * & + this % scatter % energy(gin) % data(gout) + end if + else + xs = this % scatter % scattxs(gin) + end if + case('scatter/mult') + if (present(gout)) then + if (gout < this % scatter % gmin(gin) .or. & + gout > this % scatter % gmax(gin)) then + xs = ZERO + else + xs = this % scatter % scattxs(gin) * & + this % scatter % energy(gin) % data(gout) / & + this % scatter % mult(gin) % data(gout) + end if + else + xs = this % scatter % scattxs(gin) / & + (dot_product(this % scatter % mult(gin) % data, & + this % scatter % energy(gin) % data)) + end if + case('scatter*f_mu/mult','scatter*f_mu') + if (present(gout)) then + if (gout < this % scatter % gmin(gin) .or. & + gout > this % scatter % gmax(gin)) then + xs = ZERO + else + xs = this % scatter % scattxs(gin) * & + this % scatter % energy(gin) % data(gout) * & + this % scatter % calc_f(gin, gout, mu) + if (xstype == 'scatter*f_mu/mult') then + xs = xs / this % scatter % mult(gin) % data(gout) + end if + end if + else + xs = ZERO + ! TODO (Not likely needed) + ! (asking for f_mu without asking for a group or mu would mean the + ! user of this code wants the complete 1-outgoing group distribution + ! which Im not sure what they would do with that. + end if + case default + xs = ZERO + end select + + end function mgxsiso_get_xs + + pure function mgxsang_get_xs(this, xstype, gin, gout, uvw, mu) result(xs) + class(MgxsAngle), intent(in) :: this ! The Mgxs to initialize + character(*) , intent(in) :: xstype ! Type of xs requested + integer, intent(in) :: gin ! Incoming Energy group + integer, optional, intent(in) :: gout ! Outgoing Energy group + real(8), optional, intent(in) :: uvw(3) ! Requested Angle + real(8), optional, intent(in) :: mu ! Change in angle + real(8) :: xs ! Requested x/s + + integer :: iazi, ipol + + if (present(uvw)) then + call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + select case(xstype) + case('total') + xs = this % total(gin, iazi, ipol) + case('absorption') + xs = this % absorption(gin, iazi, ipol) + case('fission') + if (allocated(this % fission)) then + xs = this % fission(gin, iazi, ipol) + else + xs = ZERO + end if + case('kappa_fission') + if (allocated(this % k_fission)) then + xs = this % k_fission(gin, iazi, ipol) + else + xs = ZERO + end if + case('nu_fission') + xs = this % nu_fission(gin, iazi, ipol) + case('chi') + if (present(gout)) then + xs = this % chi(gout, gin, iazi, ipol) + else + ! 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') + if (present(gout)) then + if (gout < this % scatter(iazi, ipol) % obj % gmin(gin) .or. & + gout > this % scatter(iazi, ipol) % obj % gmax(gin)) then + xs = ZERO + else + xs = this % scatter(iazi, ipol) % obj % scattxs(gin) * & + this % scatter(iazi, ipol) % obj % energy(gin) % data(gout) + end if + else + xs = this % scatter(iazi, ipol) % obj % scattxs(gin) + end if + case('scatter/mult') + if (present(gout)) then + if (gout < this % scatter(iazi, ipol) % obj % gmin(gin) .or. & + gout > this % scatter(iazi, ipol) % obj % gmax(gin)) then + xs = ZERO + else + xs = this % scatter(iazi, ipol) % obj % scattxs(gin) * & + this % scatter(iazi, ipol) % obj % energy(gin) % data(gout) / & + this % scatter(iazi, ipol) % obj % mult(gin) % data(gout) + end if + else + xs = this % scatter(iazi, ipol) % obj % scattxs(gin) / & + (dot_product(this % scatter(iazi, ipol) % obj % mult(gin) % data, & + this % scatter(iazi, ipol) % obj % energy(gin) % data)) + end if + case('scatter*f_mu/mult','scatter*f_mu') + if (present(gout)) then + if (gout < this % scatter(iazi, ipol) % obj % gmin(gin) .or. & + gout > this % scatter(iazi, ipol) % obj % gmax(gin)) then + xs = ZERO + else + xs = this % scatter(iazi, ipol) % obj % scattxs(gin) * & + this % scatter(iazi, ipol) % obj % energy(gin) % data(gout) + xs = xs * this % scatter(iazi, ipol) % obj % calc_f(gin, gout, mu) + if (xstype == 'scatter*f_mu/mult') then + xs = xs / & + this % scatter(iazi, ipol) % obj % mult(gin) % data(gout) + end if + end if + else + xs = ZERO + ! TODO (Not likely needed) + ! (asking for f_mu without asking for a group or mu would mean the + ! user of this code wants the complete 1-outgoing group distribution + ! which Im not sure what they would do with that. + end if + case default + xs = ZERO + end select + else + xs = ZERO + end if + + end function mgxsang_get_xs + +!=============================================================================== +! MACROXS*_COMBINE Builds a macroscopic Mgxs object from microscopic Mgxs +! objects +!=============================================================================== + + 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 + integer, intent(in) :: i_listing ! Index in listings + + ! Fill in meta-data from material information + if (mat % name == "") then + this % name = trim(to_str(mat % id)) + else + this % name = mat % name + end if + this % zaid = -mat % id + this % listing = i_listing + this % fissionable = mat % fissionable + this % scatt_type = scatt_type + + ! The following info we should initialize, but we dont need it nor + ! does it have guaranteed meaning. + this % awr = -ONE + this % kT = -ONE + + end subroutine mgxs_combine + + 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 + type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from + integer, intent(in) :: groups ! Number of E groups + integer, intent(in) :: max_order ! Maximum requested order + integer, intent(in) :: scatt_type ! How is data presented + integer, intent(in) :: i_listing ! Index in listings + + integer :: i ! loop index over nuclides + integer :: gin, gout ! group indices + real(8) :: atom_density ! atom density of a nuclide + real(8) :: norm, nuscatt + integer :: mat_max_order, order, order_dim, nuc_order_dim + real(8), allocatable :: temp_mult(:, :), mult_num(:, :), mult_denom(:, :) + real(8), allocatable :: scatt_coeffs(:, :, :) + + ! Set the meta-data + call mgxs_combine(this, mat, scatt_type, i_listing) + + ! Determine the scattering type of our data and ensure all scattering orders + ! are the same. + select type(nuc => nuclides(mat % nuclide(1)) % obj) + type is (MgxsIso) + order = size(nuc % scatter % dist(1) % data, dim=1) + end select + ! If we have tabular only data, then make sure all datasets have same size + if (scatt_type == ANGLE_HISTOGRAM) then + ! Check all scattering data to ensure it is the same size + 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!") + end select + end do + ! Ok, got our order, store the dimensionality + order_dim = order + + ! Set our Scatter Object Type + allocate(ScattDataHistogram :: this % scatter) + + else if (scatt_type == ANGLE_TABULAR) then + ! Check all scattering data to ensure it is the same size + 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 tabular scattering entries must be& + & same length!") + end select + end do + ! Ok, got our order, store the dimensionality + order_dim = order + + ! Set our Scatter Object Type + allocate(ScattDataTabular :: this % scatter) + + else if (scatt_type == ANGLE_LEGENDRE) then + ! Need to determine the maximum scattering order of all data in this material + mat_max_order = 0 + 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) & + mat_max_order = size(nuc % scatter % dist(1) % data, dim=1) + end select + end do + + ! Now need to compare this material maximum scattering order with + ! the problem wide max scatt order and use whichever is lower + order = min(mat_max_order, max_order) + ! Ok, got our order, store the dimensionality + order_dim = order + 1 + + ! Set our Scatter Object Type + allocate(ScattDataLegendre :: this % scatter) + end if + + ! Allocate and initialize data needed for macro_xs(i_mat) object + allocate(this % total(groups)) + this % total(:) = ZERO + allocate(this % absorption(groups)) + this % absorption(:) = ZERO + allocate(this % fission(groups)) + this % fission(:) = ZERO + allocate(this % k_fission(groups)) + this % k_fission(:) = ZERO + allocate(this % nu_fission(groups)) + this % nu_fission(:) = ZERO + allocate(this % chi(groups,groups)) + this % chi(:, :) = ZERO + allocate(temp_mult(groups,groups)) + temp_mult(:, :) = ZERO + allocate(mult_num(groups,groups)) + mult_num(:, :) = ZERO + allocate(mult_denom(groups,groups)) + mult_denom(:, :) = ZERO + allocate(scatt_coeffs(order_dim,groups,groups)) + scatt_coeffs(:, :, :) = ZERO + + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! Perform our operations which depend upon the type + 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(:) + if (nuc % fissionable) then + 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(:) + end if + if (allocated(nuc % k_fission)) then + this % k_fission(:) = this % k_fission(:) + atom_density * & + nuc % k_fission(:) + end if + end if + + ! Get the multiplication 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'})) + ! Developed as follows: + ! mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'} + ! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'}) + ! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / + ! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'})) + ! nuscatt_{i,g,g'} can be reconstructed from scatter % energy and + ! scatter % scattxs + do gin = 1, groups + do gout = nuc % scatter % gmin(gin), nuc % scatter % gmax(gin) + nuscatt = nuc % scatter % scattxs(gin) * & + nuc % scatter % energy(gin) % data(gout) + mult_num(gout, gin) = mult_num(gout, gin) + atom_density * & + nuscatt + mult_denom(gout, gin) = mult_denom(gout,gin) + atom_density * & + nuscatt / nuc % scatter % mult(gin) % data(gout) + end do + end do + + ! Get the complete scattering matrix + nuc_order_dim = size(nuc % 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 % scatter % get_matrix(min(nuc_order_dim, order_dim)) + + type is (MgxsAngle) + call fatal_error("Invalid passing of MgxsAngle to MgxsIso object") + end select + end do + + ! Obtain temp_mult + do gin = 1, groups + do gout = 1, groups + if (mult_denom(gout, gin) > ZERO) then + temp_mult(gout, gin) = mult_num(gout, gin) / mult_denom(gout, gin) + else + temp_mult(gout, gin) = ONE + end if + end do + end do + + ! Initialize the ScattData Object + call this % scatter % init(temp_mult, scatt_coeffs) + + ! Now normalize chi + if (mat % fissionable) then + do gin = 1, groups + norm = sum(this % chi(:, gin)) + if (norm > ZERO) then + this % chi(:, gin) = this % chi(:, gin) / norm + end if + end do + end if + + ! Deallocate temporaries + deallocate(scatt_coeffs, temp_mult, mult_num, mult_denom) + + end subroutine mgxsiso_combine + + 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 + type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from + integer, intent(in) :: groups ! Number of E groups + integer, intent(in) :: max_order ! Maximum requested order + integer, intent(in) :: scatt_type ! Legendre or Tabular Scatt? + integer, intent(in) :: i_listing ! Index in listings + + integer :: i ! loop index over nuclides + integer :: gin, gout ! group indices + real(8) :: atom_density ! atom density of a nuclide + integer :: ipol, iazi, n_pol, n_azi + real(8) :: norm, nuscatt + integer :: mat_max_order, order, order_dim, nuc_order_dim + real(8), allocatable :: temp_mult(:, :, :, :), mult_num(:, :, :, :) + real(8), allocatable :: mult_denom(:, :, :, :), scatt_coeffs(:, :, :, :, :) + + ! Set the meta-data + 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 + n_pol = -1 + n_azi = -1 + do i = 1, mat % n_nuclides + select type(nuc => nuclides(mat % nuclide(i)) % obj) + type is (MgxsAngle) + if (n_pol == -1) then + n_pol = nuc % n_pol + n_azi = nuc % n_azi + allocate(this % polar(n_pol)) + this % polar(:) = nuc % polar(:) + allocate(this % azimuthal(n_azi)) + 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!") + end if + end if + end select + end do + + ! Determine the scattering type of our data and ensure all scattering orders + ! are the same. + select type(nuc => nuclides(mat % nuclide(1)) % obj) + type is (MgxsAngle) + order = size(nuc % scatter(1,1) % obj % dist(1) % data, dim=1) + end select + ! If we have tabular only data, then make sure all datasets have same size + if (scatt_type == ANGLE_HISTOGRAM) then + ! Check all scattering data to ensure it is the same size + ! order = size(nuclides(mat % nuclide(1)) % obj % scatter % data,dim=1) + 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 histogram scattering entries must be& + & same length!") + end select + end do + ! Ok, got our order, store the dimensionality + order_dim = order + + ! Set our Scatter Object Type + allocate(this % scatter(n_azi, n_pol)) + do ipol = 1, n_pol + do iazi = 1, n_azi + allocate(ScattDataHistogram :: this % scatter(iazi, ipol) % obj) + end do + end do + + else if (scatt_type == ANGLE_TABULAR) then + ! Check all scattering data to ensure it is the same size + 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!") + end select + end do + ! Ok, got our order, store the dimensionality + order_dim = order + + ! Set our Scatter Object Type + allocate(this % scatter(n_azi, n_pol)) + do ipol = 1, n_pol + do iazi = 1, n_azi + allocate(ScattDataTabular :: this % scatter(iazi, ipol) % obj) + end do + end do + + else if (scatt_type == ANGLE_LEGENDRE) then + ! Need to determine the maximum scattering order of all data in this material + mat_max_order = 0 + do i = 1, mat % n_nuclides + select type(nuc => nuclides(mat % nuclide(i)) % obj) + type is (MgxsAngle) + if (size(nuc % scatter(1,1) % obj % dist(1) % data, dim=1) > mat_max_order) & + mat_max_order = size(nuc % scatter(1,1) % obj% dist(1) % data, dim=1) + end select + end do + + ! Now need to compare this material maximum scattering order with + ! the problem wide max scatt order and use whichever is lower + order = min(mat_max_order, max_order) + ! Ok, got our order, store the dimensionality + order_dim = order + 1 + + ! Set our Scatter Object Type + allocate(this % scatter(n_azi, n_pol)) + do ipol = 1, n_pol + do iazi = 1, n_azi + allocate(ScattDataLegendre :: this % scatter(iazi, ipol) % obj) + end do + end do + end if + + ! Allocate and initialize data within macro_xs(i_mat) object + allocate(this % total(groups, n_azi, n_pol)) + this % total(:, :, :) = ZERO + allocate(this % absorption(groups, n_azi, n_pol)) + this % absorption(:, :, :) = ZERO + allocate(this % fission(groups, n_azi, n_pol)) + this % fission(:, :, :) = ZERO + allocate(this % k_fission(groups, n_azi, n_pol)) + this % k_fission(:, :, :) = ZERO + allocate(this % nu_fission(groups, n_azi, n_pol)) + this % nu_fission(:, :, :) = ZERO + allocate(this % chi(groups, groups, n_azi, n_pol)) + this % chi(:, :, :, :) = ZERO + allocate(temp_mult(groups, groups, n_azi, n_pol)) + temp_mult(:, :, :, :) = ZERO + allocate(mult_num(groups, groups, n_azi, n_pol)) + mult_num(:, :, :, :) = ZERO + allocate(mult_denom(groups, groups, n_azi, n_pol)) + mult_denom(:, :, :, :) = ZERO + allocate(scatt_coeffs(order_dim, groups, groups, n_azi, n_pol)) + scatt_coeffs(:, :, :, :, :) = ZERO + + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! 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") + 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(:, :, :) + if (nuc % fissionable) then + 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(:, :, :) + end if + if (allocated(nuc % k_fission)) then + this % k_fission(:, :, :) = this % k_fission(:, :, :) + & + atom_density * nuc % k_fission(:, :, :) + end if + end if + + ! Get the multiplication 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'})) + ! Developed as follows: + ! mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'} + ! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'}) + ! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / + ! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'})) + ! nuscatt_{i,g,g'} can be reconstructed from scatter % energy and + ! scatter % scattxs + do ipol = 1, n_pol + do iazi = 1, n_azi + do gin = 1, groups + do gout = nuc % scatter(iazi, ipol) % obj % gmin(gin), & + nuc % scatter(iazi, ipol) % obj % gmax(gin) + nuscatt = nuc % scatter(iazi, ipol) % obj % scattxs(gin) * & + nuc % scatter(iazi, ipol) % obj % energy(gin) % data(gout) + mult_num(gout, gin, iazi, ipol) = mult_num(gout, gin, iazi, ipol) + & + atom_density * nuscatt + mult_denom(gout, gin, iazi, ipol) = & + mult_denom(gout, gin, iazi, ipol) + & + atom_density * nuscatt / & + nuc % scatter(iazi, ipol) % obj % mult(gin) % data(gout) + end do + end do + end do + end do + + ! Get the complete scattering matrix + 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) = & + scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) + & + atom_density * nuc % scatter(iazi, ipol) % obj % get_matrix(& + min(nuc_order_dim, order_dim)) + end do + end do + end select + end do + + ! Obtain temp_mult + do ipol = 1, n_pol + do iazi = 1, n_azi + do gin = 1, groups + do gout = 1, groups + if (mult_denom(gout, gin, iazi, ipol) > ZERO) then + temp_mult(gout, gin, iazi, ipol) = & + mult_num(gout, gin, iazi, ipol) / & + mult_denom(gout, gin, iazi, ipol) + else + temp_mult(gout, gin, iazi, ipol) = ONE + 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 % scatter(iazi, ipol) % obj % init( & + temp_mult(:, :, iazi, ipol), scatt_coeffs(:, :, :, iazi, ipol)) + end do + end do + + ! Now normalize chi + if (mat % fissionable) then + do ipol = 1, n_pol + do iazi = 1, n_azi + do gin = 1, groups + norm = sum(this % chi(:, gin, iazi, ipol)) + if (norm > ZERO) then + this % chi(:, gin, iazi, ipol) = this % chi(:, gin, iazi, ipol) / norm + end if + end do + end do + end do + end if + + ! Deallocate temporaries for the next material + deallocate(scatt_coeffs, temp_mult) + + end subroutine mgxsang_combine + +!=============================================================================== +! MGXS*_SAMPLE_FISSION_ENERGY samples the outgoing energy from a fission event +!=============================================================================== + + function mgxsiso_sample_fission_energy(this, gin, uvw) result(gout) + class(MgxsIso), intent(in) :: this ! Data to work with + integer, intent(in) :: gin ! Incoming energy group + real(8), intent(in) :: uvw(3) ! Particle Direction + integer :: gout ! Sampled outgoing group + real(8) :: xi ! Our random number + real(8) :: prob ! Running probability + + xi = prn() + gout = 1 + prob = this % chi(gout,gin) + + do while (prob < xi) + gout = gout + 1 + prob = prob + this % chi(gout,gin) + end do + + end function mgxsiso_sample_fission_energy + + function mgxsang_sample_fission_energy(this, gin, uvw) result(gout) + class(MgxsAngle), intent(in) :: this ! Data to work with + integer, intent(in) :: gin ! Incoming energy group + real(8), intent(in) :: uvw(3) ! Particle Direction + integer :: gout ! Sampled outgoing group + real(8) :: xi ! Our random number + real(8) :: prob ! Running probability + integer :: iazi, ipol + + call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + + xi = prn() + gout = 1 + prob = this % chi(gout, gin, iazi, ipol) + + do while (prob < xi) + gout = gout + 1 + prob = prob + this % chi(gout, gin, iazi, ipol) + end do + + end function mgxsang_sample_fission_energy + +!=============================================================================== +! MGXS*_SAMPLE_SCATTER Selects outgoing energy and angle after a scatter event +!=============================================================================== + + subroutine mgxsiso_sample_scatter(this, uvw, gin, gout, mu, wgt) + class(MgxsIso), intent(in) :: this + real(8), intent(in) :: uvw(3) ! Incoming neutron direction + integer, intent(in) :: gin ! Incoming neutron group + integer, intent(out) :: gout ! Sampled outgoin group + real(8), intent(out) :: mu ! Sampled change in angle + real(8), intent(inout) :: wgt ! Particle weight + + call this % scatter % sample(gin, gout, mu, wgt) + + end subroutine mgxsiso_sample_scatter + + subroutine mgxsang_sample_scatter(this, uvw, gin, gout, mu, wgt) + class(MgxsAngle), intent(in) :: this + real(8), intent(in) :: uvw(3) ! Incoming neutron direction + integer, intent(in) :: gin ! Incoming neutron group + integer, intent(out) :: gout ! Sampled outgoin group + real(8), intent(out) :: mu ! Sampled change in angle + real(8), intent(inout) :: wgt ! Particle weight + + integer :: iazi, ipol ! Angular indices + + call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + call this % scatter(iazi, ipol) % obj % sample(gin, gout, mu, wgt) + + end subroutine mgxsang_sample_scatter + +!=============================================================================== +! MGXS*_CALCULATE_XS determines the multi-group cross sections +! for the material the particle is currently traveling through. +!=============================================================================== + + subroutine mgxsiso_calculate_xs(this, gin, uvw, xs) + class(MgxsIso), intent(in) :: this + integer, intent(in) :: gin ! Incoming neutron group + real(8), intent(in) :: uvw(3) ! Incoming neutron direction + type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data + + xs % total = this % total(gin) + xs % elastic = this % scatter % scattxs(gin) + xs % absorption = this % absorption(gin) + xs % fission = this % fission(gin) + xs % nu_fission = this % nu_fission(gin) + + end subroutine mgxsiso_calculate_xs + + subroutine mgxsang_calculate_xs(this, gin, uvw, xs) + class(MgxsAngle), intent(in) :: this + integer, intent(in) :: gin ! Incoming neutron group + real(8), intent(in) :: uvw(3) ! Incoming neutron direction + type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data + + integer :: iazi, ipol + + call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + xs % total = this % total(gin, iazi, ipol) + xs % elastic = this % scatter(iazi, ipol) % obj % scattxs(gin) + xs % absorption = this % absorption(gin, iazi, ipol) + xs % fission = this % fission(gin, iazi, ipol) + xs % nu_fission = this % nu_fission(gin, iazi, ipol) + + end subroutine mgxsang_calculate_xs + +!=============================================================================== +! find_angle finds the closest angle on the data grid and returns that index +!=============================================================================== + + pure subroutine find_angle(polar, azimuthal, uvw, i_azi, i_pol) + real(8), intent(in) :: polar(:) ! Polar angles [0,pi] + real(8), intent(in) :: azimuthal(:) ! Azi. angles [-pi,pi] + real(8), intent(in) :: uvw(3) ! Direction of motion + integer, intent(inout) :: i_pol ! Closest polar bin + integer, intent(inout) :: i_azi ! Closest azi bin + + real(8) :: my_pol, my_azi, dangle + + ! Convert uvw to polar and azi + + my_pol = acos(uvw(3)) + my_azi = atan2(uvw(2), uvw(1)) + + ! Search for equi-binned angles + dangle = PI / real(size(polar),8) + i_pol = floor(my_pol / dangle + ONE) + dangle = TWO * PI / real(size(azimuthal),8) + i_azi = floor((my_azi + PI) / dangle + ONE) + + end subroutine find_angle + +end module mgxs_header \ No newline at end of file diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 6634bcc943..92e2d7ba33 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -8,7 +8,7 @@ module nuclide_header use endf_header, only: Function1D use error, only: fatal_error, warning use list_header, only: ListInt - use math, only: evaluate_legendre, find_angle + use math, only: evaluate_legendre use product_header, only: AngleEnergyContainer use reaction_header, only: Reaction use stl_vector, only: VectorInt @@ -19,12 +19,12 @@ module nuclide_header implicit none !=============================================================================== -! Nuclide contains the base nuclidic data for a nuclide, which does not depend -! upon how the nuclear data is represented (i.e., CE, or any variant of MG). -! The extended types, NuclideCE and NuclideMG deal with the rest +! Nuclide contains the base nuclidic data for a nuclide described as needed +! for continuous-energy neutron transport. !=============================================================================== - type, abstract :: Nuclide + type :: Nuclide + ! Nuclide meta-data character(12) :: name ! name of nuclide, e.g. 92235.03c integer :: zaid ! Z and A identifier, e.g. 92235 real(8) :: awr ! Atomic Weight Ratio @@ -32,21 +32,8 @@ module nuclide_header real(8) :: kT ! temperature in MeV (k*T) ! Fission information - logical :: fissionable ! nuclide is fissionable? + logical :: fissionable ! nuclide is fissionable? - contains - procedure(nuclide_print_), deferred :: print ! Writes nuclide info - end type Nuclide - - abstract interface - subroutine nuclide_print_(this, unit) - import Nuclide - class(Nuclide),intent(in) :: this - integer, optional, intent(in) :: unit - end subroutine nuclide_print_ - end interface - - type, extends(Nuclide) :: NuclideCE ! Energy grid information integer :: n_grid ! # of nuclide grid points integer, allocatable :: grid_index(:) ! log grid mapping indices @@ -90,129 +77,14 @@ module nuclide_header ! array; used at tally-time contains - procedure :: clear => nuclidece_clear - procedure :: print => nuclidece_print - procedure :: nu => nuclidece_nu - end type NuclideCE - - type, abstract, extends(Nuclide) :: NuclideMG - ! Scattering Order Information - integer :: order ! Order of data (Scattering for NuclideIso, - ! Number of angles for all in NuclideAngle) - integer :: scatt_type ! either legendre, histogram, or tabular. - integer :: legendre_mu_points ! Number of tabular points to use to represent - ! Legendre distribs, -1 if sample with the - ! Legendres themselves - contains - procedure(nuclidemg_init_), deferred :: init ! Initialize the data - procedure(nuclidemg_get_xs_), deferred :: get_xs ! Get the requested xs - procedure(nuclidemg_calc_f_), deferred :: calc_f ! Calculates f, given mu - end type NuclideMG - - abstract interface - - subroutine nuclidemg_init_(this, node_xsdata, groups, get_kfiss, get_fiss) - import NuclideMG, Node - class(NuclideMG), intent(inout) :: this ! Working Object - type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml - integer, intent(in) :: groups ! Number of Energy groups - logical, intent(in) :: get_kfiss ! Need Kappa-Fission? - logical, intent(in) :: get_fiss ! Should we get fiss data? - end subroutine nuclidemg_init_ - - function nuclidemg_get_xs_(this, g, xstype, gout, uvw, mu, i_azi, i_pol) & - result(xs) - import NuclideMG - class(NuclideMG), intent(in) :: this - integer, intent(in) :: g ! Incoming Energy group - character(*), intent(in) :: xstype ! Cross Section Type - integer, optional, intent(in) :: gout ! Outgoing Group - real(8), optional, intent(in) :: uvw(3) ! Requested Angle - real(8), optional, intent(in) :: mu ! Change in angle - integer, optional, intent(in) :: i_azi ! Azimuthal Index - integer, optional, intent(in) :: i_pol ! Polar Index - real(8) :: xs ! Resultant xs - end function nuclidemg_get_xs_ - - pure function nuclidemg_calc_f_(this, gin, gout, mu, uvw, i_azi, i_pol) result(f) - import NuclideMG - class(NuclideMG), intent(in) :: this - integer, intent(in) :: gin ! Incoming Energy Group - integer, intent(in) :: gout ! Outgoing Energy Group - real(8), intent(in) :: mu ! Angle of interest - real(8), intent(in), optional :: uvw(3) ! Direction vector - integer, intent(in), optional :: i_azi ! Incoming Energy Group - integer, intent(in), optional :: i_pol ! Outgoing Energy Group - real(8) :: f ! Return value of f(mu) - - end function nuclidemg_calc_f_ - end interface - -!=============================================================================== -! NuclideIso contains the base MGXS data for a nuclide specifically for -! isotropically weighted MGXS -!=============================================================================== - - type, extends(NuclideMG) :: NuclideIso - - ! Microscopic cross sections - real(8), allocatable :: total(:) ! total cross section - real(8), allocatable :: absorption(:) ! absorption cross section - real(8), allocatable :: scatter(:,:,:) ! scattering information - real(8), allocatable :: nu_fission(:,:) ! fission matrix (Gout x Gin) - real(8), allocatable :: k_fission(:) ! kappa-fission - real(8), allocatable :: fission(:) ! neutron production - real(8), allocatable :: chi(:) ! Fission Spectra - real(8), allocatable :: mult(:,:) ! Scatter multiplicity (Gout x Gin) - - contains - procedure :: init => nuclideiso_init ! Initialize Nuclidic MGXS Data - procedure :: print => nuclideiso_print ! Writes nuclide info - procedure :: get_xs => nuclideiso_get_xs ! Gets Size of Data w/in Object - procedure :: calc_f => nuclideiso_calc_f ! Calcs f given mu - end type NuclideIso - -!=============================================================================== -! NuclideAngle contains the base MGXS data for a nuclide specifically for -! explicit angle-dependent weighted MGXS -!=============================================================================== - - type, extends(NuclideMG) :: NuclideAngle - - ! Microscopic cross sections. Dimensions are: (n_pol, n_azi, Nl, Ng, Ng) - real(8), allocatable :: total(:,:,:) ! total cross section - real(8), allocatable :: absorption(:,:,:) ! absorption cross section - real(8), allocatable :: scatter(:,:,:,:,:) ! scattering information - real(8), allocatable :: nu_fission(:,:,:,:) ! fission matrix (Gout x Gin) - real(8), allocatable :: k_fission(:,:,:) ! kappa-fission - real(8), allocatable :: fission(:,:,:) ! neutron production - real(8), allocatable :: chi(:,:,:) ! Fission Spectra - real(8), allocatable :: mult(:,:,:,:) ! Scatter multiplicity (Gout x Gin) - - ! In all cases, right-most indices are theta, phi - integer :: n_pol ! Number of polar angles - integer :: n_azi ! Number of azimuthal angles - real(8), allocatable :: polar(:) ! polar angles - real(8), allocatable :: azimuthal(:) ! azimuthal angles - - contains - procedure :: init => nuclideangle_init ! Initialize Nuclidic MGXS Data - procedure :: print => nuclideangle_print ! Gets Size of Data w/in Object - procedure :: get_xs => nuclideangle_get_xs ! Gets Size of Data w/in Object - procedure :: calc_f => nuclideangle_calc_f ! Calcs f given mu - end type NuclideAngle - -!=============================================================================== -! NUCLIDEMGCONTAINER pointer array for storing Nuclides -!=============================================================================== - - type NuclideMGContainer - class(NuclideMG), pointer :: obj - end type NuclideMGContainer + procedure :: clear => nuclide_clear + procedure :: print => nuclide_print + procedure :: nu => nuclide_nu + end type Nuclide !=============================================================================== ! NUCLIDE0K temporarily contains all 0K cross section data and other parameters -! needed to treat resonance scattering before transferring them to NuclideCE +! needed to treat resonance scattering before transferring them to Nuclide !=============================================================================== type Nuclide0K @@ -285,406 +157,27 @@ module nuclide_header contains !=============================================================================== -! NUCLIDE_*_INIT reads in the data from the XML file, as already accessed +! NUCLIDE_CLEAR resets and deallocates data in Nuclide !=============================================================================== - subroutine nuclidemg_init(this, node_xsdata) - class(NuclideMG), intent(inout) :: this ! Working Object - type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml + subroutine nuclide_clear(this) - type(Node), pointer :: node_legendre_mu - character(MAX_LINE_LEN) :: temp_str - logical :: enable_leg_mu + class(Nuclide), intent(inout) :: this ! The Nuclide object to clear - ! Load the data - call get_node_value(node_xsdata, "name", this % name) - this % name = to_lower(this % name) - if (check_for_node(node_xsdata, "kT")) then - call get_node_value(node_xsdata, "kT", this % kT) - else - this % kT = ZERO - end if - if (check_for_node(node_xsdata, "zaid")) then - call get_node_value(node_xsdata, "zaid", this % zaid) - else - this % zaid = -1 - end if - if (check_for_node(node_xsdata, "scatt_type")) then - call get_node_value(node_xsdata, "scatt_type", temp_str) - temp_str = trim(to_lower(temp_str)) - if (temp_str == 'legendre') then - this % scatt_type = ANGLE_LEGENDRE - else if (temp_str == 'histogram') then - this % scatt_type = ANGLE_HISTOGRAM - else if (temp_str == 'tabular') then - this % scatt_type = ANGLE_TABULAR - else - call fatal_error("Invalid Scatt Type Option!") - end if - else - this % scatt_type = ANGLE_LEGENDRE - end if + integer :: i ! Loop counter - if (check_for_node(node_xsdata, "order")) then - call get_node_value(node_xsdata, "order", this % order) - else - call fatal_error("Order Must Be Provided!") - end if + if (associated(this % urr_data)) deallocate(this % urr_data) - ! Get scattering treatment - if (check_for_node(node_xsdata, "tabular_legendre")) then - call get_node_ptr(node_xsdata, "tabular_legendre", node_legendre_mu) - if (check_for_node(node_legendre_mu, "enable")) then - call get_node_value(node_legendre_mu, "enable", temp_str) - temp_str = trim(to_lower(temp_str)) - if (temp_str == 'true' .or. temp_str == '1') then - enable_leg_mu = .true. - elseif (temp_str == 'false' .or. temp_str == '0') then - enable_leg_mu = .false. - this % legendre_mu_points = 1 - else - call fatal_error("Unrecognized tabular_legendre/enable: " // temp_str) - end if - else - enable_leg_mu = .true. - this % legendre_mu_points = 33 - end if - if (enable_leg_mu .and. & - check_for_node(node_legendre_mu, "num_points")) then - call get_node_value(node_legendre_mu, "num_points", & - this % legendre_mu_points) - if (this % legendre_mu_points <= 0) then - call fatal_error("num_points element must be positive and non-zero!") - end if - this % legendre_mu_points = -1 * this % legendre_mu_points - end if - else - this % legendre_mu_points = 1 - end if + call this % reaction_index % clear() - if (check_for_node(node_xsdata, "fissionable")) then - call get_node_value(node_xsdata, "fissionable", temp_str) - temp_str = to_lower(temp_str) - if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then - this % fissionable = .true. - else - this % fissionable = .false. - end if - else - call fatal_error("Fissionable element must be set!") - end if - - end subroutine nuclidemg_init - - subroutine nuclideiso_init(this, node_xsdata, groups, get_kfiss, get_fiss) - class(NuclideIso), intent(inout) :: this ! Working Object - type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml - integer, intent(in) :: groups ! Number of Energy groups - logical, intent(in) :: get_kfiss ! Need Kappa-Fission? - logical, intent(in) :: get_fiss ! Need fiss data? - - real(8), allocatable :: temp_arr(:) - integer :: arr_len - integer :: order_dim - - ! Call generic data gathering routine - call nuclidemg_init(this, node_xsdata) - - ! Load the more specific data - if (this % fissionable) then - - if (check_for_node(node_xsdata, "chi")) then - ! Get chi - allocate(this % chi(groups)) - call get_node_array(node_xsdata, "chi", this % chi) - - ! Get nu_fission (as a vector) - if (check_for_node(node_xsdata, "nu_fission")) then - allocate(temp_arr(groups * 1)) - call get_node_array(node_xsdata, "nu_fission", temp_arr) - allocate(this % nu_fission(groups, 1)) - this % nu_fission = reshape(temp_arr, (/groups, 1/)) - deallocate(temp_arr) - else - call fatal_error("If fissionable, must provide nu_fission!") - end if - - else - ! Get nu_fission (as a matrix) - if (check_for_node(node_xsdata, "nu_fission")) then - - allocate(temp_arr(groups*groups)) - call get_node_array(node_xsdata, "nu_fission", temp_arr) - allocate(this % nu_fission(groups, groups)) - this % nu_fission = reshape(temp_arr, (/groups, groups/)) - deallocate(temp_arr) - else - call fatal_error("If fissionable, must provide nu_fission!") - end if - end if - if (get_fiss) then - allocate(this % fission(groups)) - if (check_for_node(node_xsdata, "fission")) then - call get_node_array(node_xsdata, "fission", this % fission) - else - call fatal_error("Fission data missing, required due to fission& - & tallies in tallies.xml file!") - end if - end if - if (get_kfiss) then - allocate(this % k_fission(groups)) - if (check_for_node(node_xsdata, "kappa_fission")) then - call get_node_array(node_xsdata, "kappa_fission", this % k_fission) - else - call fatal_error("kappa_fission data missing, required due to & - &kappa-fission tallies in tallies.xml file!") - end if - end if - end if - - allocate(this % absorption(groups)) - if (check_for_node(node_xsdata, "absorption")) then - call get_node_array(node_xsdata, "absorption", this % absorption) - else - call fatal_error("Must provide absorption!") - end if - - if (this % scatt_type == ANGLE_LEGENDRE) then - order_dim = this % order + 1 - else if (this % scatt_type == ANGLE_HISTOGRAM) then - order_dim = this % order - else if (this % scatt_type == ANGLE_TABULAR) then - order_dim = this % order - end if - - allocate(this % scatter(groups, groups, order_dim)) - if (check_for_node(node_xsdata, "scatter")) then - allocate(temp_arr(groups * groups * order_dim)) - call get_node_array(node_xsdata, "scatter", temp_arr) - this % scatter = reshape(temp_arr, (/groups, groups, order_dim/)) - deallocate(temp_arr) - else - call fatal_error("Must provide scatter!") - return - end if - - - allocate(this % total(groups)) - if (check_for_node(node_xsdata, "total")) then - call get_node_array(node_xsdata, "total", this % total) - else - this % total = this % absorption + sum(this%scatter(:,:,1),dim=1) - end if - - ! Get Mult Data - allocate(this % mult(groups, groups)) - if (check_for_node(node_xsdata, "multiplicity")) then - arr_len = get_arraysize_double(node_xsdata, "multiplicity") - if (arr_len == groups * groups) then - allocate(temp_arr(arr_len)) - call get_node_array(node_xsdata, "multiplicity", temp_arr) - this % mult = reshape(temp_arr, (/groups, groups/)) - deallocate(temp_arr) - else - call fatal_error("Multiplicity length not same as number of groups& - & squared!") - return - end if - else - this % mult = ONE - end if - - end subroutine nuclideiso_init - - subroutine nuclideangle_init(this, node_xsdata, groups, get_kfiss, get_fiss) - class(NuclideAngle), intent(inout) :: this ! Working Object - type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml - integer, intent(in) :: groups ! Number of Energy groups - logical, intent(in) :: get_kfiss ! Need Kappa-Fission? - logical, intent(in) :: get_fiss ! Should we get fiss data? - - real(8), allocatable :: temp_arr(:) - integer :: arr_len - real(8) :: dangle - integer :: iangle - integer :: order_dim - - ! Call generic data gathering routine - call nuclidemg_init(this, node_xsdata) - - if (this % scatt_type == ANGLE_LEGENDRE) then - order_dim = this % order + 1 - else if (this % scatt_type == ANGLE_HISTOGRAM) then - order_dim = this % order - else if (this % scatt_type == ANGLE_TABULAR) then - order_dim = this % order - end if - - 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!") - 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!") - end if - - ! Load angle data, if present (else equally spaced) - allocate(this % polar(this % n_pol)) - allocate(this % azimuthal(this % n_azi)) - if (check_for_node(node_xsdata, "polar")) then - call fatal_error("User-Specified polar angle bins not yet supported!") - ! When this feature is supported, this line will be activated - call get_node_array(node_xsdata, "polar", this % polar) - else - dangle = PI / real(this % n_pol,8) - do iangle = 1, this % n_pol - this % polar(iangle) = (real(iangle,8) - HALF) * dangle - end do - end if - if (check_for_node(node_xsdata, "azimuthal")) then - call fatal_error("User-Specified azimuthal angle bins not yet supported!") - ! When this feature is supported, this line will be activated - call get_node_array(node_xsdata, "azimuthal", this % azimuthal) - else - dangle = TWO * PI / real(this % n_azi,8) - do iangle = 1, this % n_azi - this % azimuthal(iangle) = -PI + (real(iangle,8) - HALF) * dangle - end do - end if - - ! Load the more specific data - if (this % fissionable) then - - if (check_for_node(node_xsdata, "chi")) then - ! Get chi - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call get_node_array(node_xsdata, "chi", temp_arr) - allocate(this % chi(groups, this % n_azi, this % n_pol)) - this % chi = reshape(temp_arr, (/groups, this % n_azi, this % n_pol/)) - deallocate(temp_arr) - - ! Get nu_fission (as a vector) - 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) - allocate(this % nu_fission(groups, 1, this % n_azi, this % n_pol)) - this % nu_fission = reshape(temp_arr, (/groups, 1, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) - else - call fatal_error("If fissionable, must provide nu_fission!") - end if - - else - ! Get nu_fission (as a matrix) - 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) - allocate(this % nu_fission(groups, groups, this % n_azi, this % n_pol)) - this % nu_fission = reshape(temp_arr, (/groups, groups, & - this % n_azi, this % n_pol/)) - deallocate(temp_arr) - else - call fatal_error("If fissionable, must provide nu_fission!") - end if - end if - if (get_fiss) then - if (check_for_node(node_xsdata, "fission")) then - 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/)) - deallocate(temp_arr) - else - call fatal_error("Fission data missing, required due to fission& - & tallies in tallies.xml file!") - end if - end if - if (get_kfiss) then - if (check_for_node(node_xsdata, "kappa_fission")) then - 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/)) - deallocate(temp_arr) - else - call fatal_error("kappa_fission data missing, required due to & - &kappa-fission tallies in tallies.xml file!") - end if - end if - 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/)) - deallocate(temp_arr) - else - call fatal_error("Must provide absorption!") - end if - - allocate(this % scatter(groups, groups, order_dim, this % n_azi, this % n_pol)) - if (check_for_node(node_xsdata, "scatter")) then - allocate(temp_arr(groups * groups * order_dim * this % n_azi * this%n_pol)) - call get_node_array(node_xsdata, "scatter", temp_arr) - this % scatter = reshape(temp_arr, (/groups, groups, order_dim, & - this%n_azi,this%n_pol/)) - deallocate(temp_arr) - else - call fatal_error("Must provide scatter!") - end if - - 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) - allocate(this % total(groups, this % n_azi, this % n_pol)) - this % total = reshape(temp_arr, (/groups, this % n_azi, this % n_pol/)) - deallocate(temp_arr) - else - this % total = this % absorption + sum(this%scatter(:,:,1,:,:),dim=1) - end if - - ! Get Mult Data - allocate(this % mult(groups, groups, this % n_azi, this % n_pol)) - if (check_for_node(node_xsdata, "multiplicity")) then - arr_len = get_arraysize_double(node_xsdata, "multiplicity") - 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) - this % mult = reshape(temp_arr, (/groups, groups, this % n_azi, this % n_pol/)) - deallocate(temp_arr) - else - call fatal_error("Multiplicity Length Does Not Match!") - end if - else - this % mult = ONE - end if - - end subroutine nuclideangle_init + end subroutine nuclide_clear !=============================================================================== -! NUCLIDECE_CLEAR resets and deallocates data in Nuclide, NuclideIso -! or NuclideAngle +! NUCLIDE_NU is an interface to the number of fission neutrons produced !=============================================================================== - subroutine nuclidece_clear(this) - - class(NuclideCE), intent(inout) :: this ! The Nuclide object to clear - - if (associated(this % urr_data)) deallocate(this % urr_data) - - call this % reaction_index % clear() - - end subroutine nuclidece_clear - - function nuclidece_nu(this, E, emission_mode, group) result(nu) - class(NuclideCE), intent(in) :: this + function nuclide_nu(this, E, emission_mode, group) result(nu) + class(Nuclide), intent(in) :: this real(8), intent(in) :: E integer, intent(in) :: emission_mode integer, optional, intent(in) :: group @@ -742,15 +235,16 @@ module nuclide_header end if end select - end function nuclidece_nu + end function nuclide_nu + !=============================================================================== ! NUCLIDE*_PRINT displays information about a continuous-energy neutron ! cross_section table and its reactions and secondary angle/energy distributions !=============================================================================== - subroutine nuclidece_print(this, unit) - class(NuclideCE), intent(in) :: this + subroutine nuclide_print(this, unit) + class(Nuclide), intent(in) :: this integer, intent(in), optional :: unit integer :: i ! loop index over nuclides @@ -823,346 +317,6 @@ module nuclide_header ! Blank line at end of nuclide write(unit_,*) - end subroutine nuclidece_print - - subroutine nuclidemg_print(this, unit_) - class(NuclideMG), intent(in) :: this - integer, intent(in) :: unit_ - - character(MAX_LINE_LEN) :: temp_str - - ! Basic nuclide information - write(unit_,*) 'Nuclide ' // trim(this % name) - if (this % zaid > 0) then - ! Dont print if data was macroscopic and thus zaid & AWR would be nonsense - write(unit_,*) ' zaid = ' // trim(to_str(this % zaid)) - write(unit_,*) ' awr = ' // trim(to_str(this % awr)) - end if - write(unit_,*) ' kT = ' // trim(to_str(this % kT)) - if (this % scatt_type == ANGLE_LEGENDRE) then - temp_str = "Legendre" - write(unit_,*) ' Scattering Type = ' // trim(temp_str) - write(unit_,*) ' # of Scatter Moments = ' // & - trim(to_str(this % order)) - else if (this % scatt_type == ANGLE_HISTOGRAM) then - temp_str = "Histogram" - write(unit_,*) ' Scattering Type = ' // trim(temp_str) - write(unit_,*) ' # of Scatter Bins = ' // & - trim(to_str(this % order)) - else if (this % scatt_type == ANGLE_TABULAR) then - temp_str = "Tabular" - write(unit_,*) ' Scattering Type = ' // trim(temp_str) - write(unit_,*) ' # of Scatter Points = ' // trim(to_str(this % order)) - end if - write(unit_,*) ' Fissionable = ', this % fissionable - - end subroutine nuclidemg_print - - subroutine nuclideiso_print(this, unit) - - class(NuclideIso), intent(in) :: this - integer, optional, intent(in) :: unit - - integer :: unit_ ! unit to write to - integer :: size_total, size_scattmat, size_mgxs - - ! set default unit for writing information - if (present(unit)) then - unit_ = unit - else - unit_ = OUTPUT_UNIT - end if - - ! Write Basic Nuclide Information - call nuclidemg_print(this, unit_) - - ! Determine size of mgxs and scattering matrices - size_scattmat = (size(this % scatter) + size(this % mult)) * 8 - size_mgxs = size(this % total) + size(this % absorption) + & - size(this % nu_fission) + size(this % k_fission) + & - size(this % fission) + size(this % chi) - size_mgxs = size_mgxs * 8 - - ! Calculate total memory - size_total = size_scattmat + size_mgxs - - ! Write memory used - write(unit_,*) ' Memory Requirements' - write(unit_,*) ' Cross sections = ' // trim(to_str(size_mgxs)) // ' bytes' - write(unit_,*) ' Scattering Matrices = ' // & - trim(to_str(size_scattmat)) // ' bytes' - write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' - - ! Blank line at end of nuclide - write(unit_,*) - - end subroutine nuclideiso_print - - subroutine nuclideangle_print(this, unit) - - class(NuclideAngle), intent(in) :: this - integer, optional, intent(in) :: unit - - integer :: unit_ ! unit to write to - integer :: size_total, size_scattmat, size_mgxs - - ! set default unit for writing information - if (present(unit)) then - unit_ = unit - else - unit_ = OUTPUT_UNIT - end if - - ! Write Basic Nuclide Information - call nuclidemg_print(this, unit_) - write(unit_,*) ' # of Polar Angles = ' // trim(to_str(this % n_pol)) - write(unit_,*) ' # of Azimuthal Angles = ' // trim(to_str(this % n_azi)) - - ! Determine size of mgxs and scattering matrices - size_scattmat = (size(this % scatter) + size(this % mult)) * 8 - size_mgxs = size(this % total) + size(this % absorption) + & - size(this % nu_fission) + size(this % k_fission) + & - size(this % fission) + size(this % chi) - size_mgxs = size_mgxs * 8 - - ! Calculate total memory - size_total = size_scattmat + size_mgxs - - ! Write memory used - write(unit_,*) ' Memory Requirements' - write(unit_,*) ' Cross sections = ' // trim(to_str(size_mgxs)) // ' bytes' - write(unit_,*) ' Scattering Matrices = ' // & - trim(to_str(size_scattmat)) // ' bytes' - write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' - - ! Blank line at end of nuclide - write(unit_,*) - - - end subroutine nuclideangle_print - -!=============================================================================== -! NUCLIDE*_GET_XS Returns the requested data type -!=============================================================================== - - function nuclideiso_get_xs(this, g, xstype, gout, uvw, mu, i_azi, i_pol) & - result(xs) - class(NuclideIso), intent(in) :: this - integer, intent(in) :: g ! Incoming Energy group - character(*), intent(in) :: xstype ! Cross Section Type - integer, optional, intent(in) :: gout ! Outgoing Group - real(8), optional, intent(in) :: uvw(3) ! Requested Angle - real(8), optional, intent(in) :: mu ! Change in angle - integer, optional, intent(in) :: i_azi ! Azimuthal Index - integer, optional, intent(in) :: i_pol ! Polar Index - real(8) :: xs ! Resultant xs - - xs = ZERO - - if ((xstype == 'nu_fission' .or. xstype == 'fission' .or. xstype =='chi' & - .or. xstype =='k_fission') .and. (.not. this % fissionable)) then - return - end if - - if (present(gout)) then - select case(xstype) - case('mult') - xs = this % mult(gout,g) - case('nu_fission') - xs = this % nu_fission(gout,g) - case('f_mu', 'f_mu/mult') - xs = this % calc_f(g, gout, mu) - if (xstype == 'f_mu/mult') then - xs = xs / this % mult(gout,g) - end if - end select - else - select case(xstype) - case('total') - xs = this % total(g) - case('absorption') - xs = this % absorption(g) - case('fission') - xs = this % fission(g) - case('k_fission') - if (allocated(this % k_fission)) then - xs = this % k_fission(g) - end if - case('chi') - xs = this % chi(g) - case('scatter') - xs = this % total(g) - this % absorption(g) - end select - end if - end function nuclideiso_get_xs - - function nuclideangle_get_xs(this, g, xstype, gout, uvw, mu, i_azi, i_pol) & - result(xs) - class(NuclideAngle), intent(in) :: this - integer, intent(in) :: g ! Incoming Energy group - character(*), intent(in) :: xstype ! Cross Section Type - integer, optional, intent(in) :: gout ! Outgoing Group - real(8), optional, intent(in) :: mu ! Change in angle - real(8), optional, intent(in) :: uvw(3) ! Requested Angle - integer, optional, intent(in) :: i_azi ! Azimuthal Index - integer, optional, intent(in) :: i_pol ! Polar Index - real(8) :: xs ! Resultant xs - - integer :: i_azi_, i_pol_ - - xs = ZERO - - if ((xstype == 'nu_fission' .or. xstype == 'fission' .or. xstype =='chi' & - .or. xstype =='k_fission') .and. (.not. this % fissionable)) then - return - end if - - if (present(i_azi) .and. present(i_pol)) then - i_azi_ = i_azi - i_pol_ = i_pol - else - call find_angle(this % polar, this % azimuthal, uvw, i_azi_, i_pol_) - end if - - if (present(gout)) then - select case(xstype) - case('mult') - xs = this % mult(gout,g,i_azi_,i_pol_) - case('nu_fission') - xs = this % nu_fission(gout,g,i_azi_,i_pol_) - case('chi') - xs = this % chi(gout,i_azi_,i_pol_) - case('f_mu', 'f_mu/mult') - xs = this % calc_f(g, gout, mu, I_AZI=i_azi_, I_POL=i_pol_) - if (xstype == 'f_mu/mult') then - xs = xs / this % mult(gout,g,i_azi_,i_pol_) - end if - end select - else - select case(xstype) - case('total') - xs = this % total(g,i_azi_,i_pol_) - case('absorption') - xs = this % absorption(g,i_azi_,i_pol_) - case('fission') - xs = this % fission(g,i_azi_,i_pol_) - case('k_fission') - if (allocated(this % k_fission)) then - xs = this % k_fission(g,i_azi_,i_pol_) - end if - case('chi') - xs = this % chi(g,i_azi_,i_pol_) - case('scatter') - xs = this % total(g,i_azi_,i_pol_) - this % absorption(g,i_azi_,i_pol_) - end select - end if - - end function nuclideangle_get_xs - -!=============================================================================== -! NUCLIDE*_CALC_F Finds the value of f(mu), the scattering angle probability, -! given mu -!=============================================================================== - - pure function nuclideiso_calc_f(this, gin, gout, mu, uvw, i_azi, i_pol) & - result(f) - class(NuclideIso), intent(in) :: this - integer, intent(in) :: gin ! Incoming Energy Group - integer, intent(in) :: gout ! Outgoing Energy Group - real(8), intent(in) :: mu ! Angle of interest - real(8), intent(in), optional :: uvw(3) ! Direction vector - integer, intent(in), optional :: i_azi ! Incoming Energy Group - integer, intent(in), optional :: i_pol ! Outgoing Energy Group - real(8) :: f ! Return value of f(mu) - - real(8) :: dmu, r - integer :: imu - - if (this % scatt_type == ANGLE_LEGENDRE) then - f = evaluate_legendre(this % scatter(gout,gin,:), mu) - else if (this % scatt_type == ANGLE_TABULAR) then - dmu = TWO / real(this % order - 1,8) - ! Find mu bin algebraically, knowing that the spacing is equal - f = (mu + ONE) / dmu + ONE - imu = floor(f) - ! But save the amount that mu is past the previous index - ! so we can use interpolation later. - f = f - real(imu,8) - ! Adjust so interpolation works on the last bin if necessary - if (imu == size(this % scatter, dim=3)) then - imu = imu - 1 - end if - - ! Now intepolate to find f(mu) - r = f / dmu - f = (ONE - r) * this % scatter(gout,gin,imu) + & - r * this % scatter(gout,gin,imu+1) - else ! (ANGLE_HISTOGRAM) - dmu = TWO / real(this % order,8) - ! Find mu bin algebraically, knowing that the spacing is equal - imu = floor((mu + ONE) / dmu + ONE) - ! Adjust so interpolation works on the last bin if necessary - if (imu == size(this % scatter, dim=3)) then - imu = imu - 1 - end if - f = this % scatter(gout, gin, imu) - - end if - - end function nuclideiso_calc_f - - pure function nuclideangle_calc_f(this, gin, gout, mu, uvw, i_azi, & - i_pol) result(f) - class(NuclideAngle), intent(in) :: this - integer, intent(in) :: gin ! Incoming Energy Group - integer, intent(in) :: gout ! Outgoing Energy Group - real(8), intent(in) :: mu ! Angle of interest - real(8), intent(in), optional :: uvw(3) ! Direction vector - integer, intent(in), optional :: i_azi ! Incoming Energy Group - integer, intent(in), optional :: i_pol ! Outgoing Energy Group - real(8) :: f ! Return value of f(mu) - - real(8) :: dmu, r - integer :: imu - integer :: i_azi_, i_pol_ - if (present(i_azi) .and. present(i_pol)) then - i_azi_ = i_azi - i_pol_ = i_pol - else if (present(uvw)) then - call find_angle(this % polar, this % azimuthal, uvw, i_azi_, i_pol_) - end if - - if (this % scatt_type == ANGLE_LEGENDRE) then - f = evaluate_legendre(this % scatter(gout,gin,:,i_azi_,i_pol_), mu) - else if (this % scatt_type == ANGLE_TABULAR) then - dmu = TWO / real(this % order - 1,8) - ! Find mu bin algebraically, knowing that the spacing is equal - f = (mu + ONE) / dmu + ONE - imu = floor(f) - ! But save the amount that mu is past the previous index - ! so we can use interpolation later. - f = f - real(imu,8) - ! Adjust so interpolation works on the last bin if necessary - if (imu == size(this % scatter, dim=3)) then - imu = imu - 1 - end if - - ! Now intepolate to find f(mu) - r = f / dmu - f = (ONE - r) * this % scatter(gout,gin,imu,i_azi_,i_pol_) + & - r * this % scatter(gout,gin,imu+1,i_azi_,i_pol_) - else ! (ANGLE_HISTOGRAM) - dmu = TWO / real(this % order,8) - ! Find mu bin algebraically, knowing that the spacing is equal - imu = floor((mu + ONE) / dmu + ONE) - ! Adjust so interpolation works on the last bin if necessary - if (imu == size(this % scatter, dim=3)) then - imu = imu - 1 - end if - f = this % scatter(gout, gin, imu,i_azi_,i_pol_) - - end if - - end function nuclideangle_calc_f + end subroutine nuclide_print end module nuclide_header diff --git a/src/output.F90 b/src/output.F90 index 786d9a10e3..d56de1f3d2 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -335,10 +335,10 @@ contains ! Open log file for writing open(NEWUNIT=unit_xs, FILE=path, STATUS='replace', ACTION='write') - ! Write header - call header("CROSS SECTION TABLES", unit=unit_xs) - if (run_CE) then + ! Write header + call header("CROSS SECTION TABLES", unit=unit_xs) + NUCLIDE_LOOP: do i = 1, n_nuclides_total ! Print information about nuclide call nuclides(i) % print(unit=unit_xs) @@ -349,10 +349,17 @@ contains call sab_tables(i) % print(unit=unit_xs) end do SAB_TABLES_LOOP else + ! Write header + call header("MGXS LIBRARY TABLES", unit=unit_xs) NuclideMG_LOOP: do i = 1, n_nuclides_total ! Print information about nuclide call nuclides_mg(i) % obj % print(unit=unit_xs) end do NuclideMG_LOOP + call header("MATERIAL MGXS TABLES", unit=unit_xs) + MATERIAL_LOOP: do i = 1, n_materials + ! Print information about Materials + call macro_xs(i) % obj % print(unit=unit_xs) + end do MATERIAL_LOOP end if ! Close cross section summary file @@ -903,7 +910,11 @@ contains write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & "Total Material" else - i_listing = nuclides(i_nuclide) % listing + if (run_CE) then + i_listing = nuclides(i_nuclide) % listing + else + i_listing = nuclides_MG(i_nuclide) % obj % listing + end if write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & trim(xs_listings(i_listing) % alias) end if diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 4ad4119b76..8544cb38bb 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -202,7 +202,7 @@ contains this % last_g = int(src % E) this % E = energy_bin_avg(this % g) end if - this % last_E = src % E + this % last_E = this % E end subroutine initialize_from_source diff --git a/src/physics.F90 b/src/physics.F90 index d6c4c45b05..09a43aa99f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -80,7 +80,7 @@ contains integer :: i_nuclide ! index in nuclides array integer :: i_nuc_mat ! index in material's nuclides array integer :: i_reaction ! index in nuc % reactions array - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat) @@ -203,7 +203,7 @@ contains real(8) :: f real(8) :: prob real(8) :: cutoff - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc ! Get pointer to nuclide nuc => nuclides(i_nuclide) @@ -301,7 +301,7 @@ contains real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc ! copy incoming direction uvw_old(:) = p % coord(1) % uvw @@ -416,7 +416,7 @@ contains real(8) :: v_cm(3) ! velocity of center-of-mass real(8) :: v_t(3) ! velocity of target nucleus real(8) :: uvw_cm(3) ! directional cosines in center-of-mass - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc ! get pointer to nuclide nuc => nuclides(i_nuclide) @@ -742,7 +742,7 @@ contains !=============================================================================== subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) - type(NuclideCE), intent(in) :: nuc ! target nuclide at temperature T + type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T real(8), intent(out) :: v_target(3) ! target velocity real(8), intent(in) :: v_neut(3) ! neutron velocity real(8), intent(in) :: E ! particle energy @@ -987,7 +987,7 @@ contains !=============================================================================== subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw) - type(NuclideCE), intent(in) :: nuc ! target nuclide at temperature + type(Nuclide), intent(in) :: nuc ! target nuclide at temperature real(8), intent(out) :: v_target(3) real(8), intent(in) :: E real(8), intent(in) :: uvw(3) @@ -1073,7 +1073,7 @@ contains real(8) :: nu_t ! total nu real(8) :: weight ! weight adjustment for ufs method logical :: in_mesh ! source site in ufs mesh? - type(NuclideCE), pointer :: nuc + type(Nuclide), pointer :: nuc ! Get pointers nuc => nuclides(i_nuclide) @@ -1169,10 +1169,10 @@ contains !=============================================================================== subroutine sample_fission_neutron(nuc, rxn, E_in, site) - type(NuclideCE), intent(in) :: nuc - type(Reaction), intent(in) :: rxn - real(8), intent(in) :: E_in - type(Bank), intent(inout) :: site + type(Nuclide), intent(in) :: nuc + type(Reaction), intent(in) :: rxn + real(8), intent(in) :: E_in + type(Bank), intent(inout) :: site integer :: group ! index on nu energy grid / precursor group integer :: n_sample ! number of resamples @@ -1278,7 +1278,7 @@ contains !=============================================================================== subroutine inelastic_scatter(nuc, rxn, p) - type(NuclideCE), intent(in) :: nuc + type(Nuclide), intent(in) :: nuc type(Reaction), intent(in) :: rxn type(Particle), intent(inout) :: p diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 6a58540c17..2e5e467c14 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -5,9 +5,9 @@ module physics_mg use constants use error, only: fatal_error, warning use global - use macroxs_header, only: MacroXS, MacroXSContainer use material_header, only: Material use math, only: rotate_angle + use mgxs_header, only: Mgxs, MgxsContainer use mesh, only: get_mesh_indices use output, only: write_message use particle_header, only: Particle @@ -77,6 +77,7 @@ contains call create_fission_sites(p, p % secondary_bank, p % n_secondary) end if end if + ! If survival biasing is being used, the following subroutine adjusts the ! weight of the particle. Otherwise, it checks to see if absorption occurs @@ -178,7 +179,7 @@ contains real(8) :: phi ! fission neutron azimuthal angle real(8) :: weight ! weight adjustment for ufs method logical :: in_mesh ! source site in ufs mesh? - class(MacroXS), pointer :: xs + class(Mgxs), pointer :: xs ! Get Pointers xs => macro_xs(p % material) % obj @@ -241,14 +242,12 @@ contains ! Set weight of fission bank site bank_array(i) % wgt = ONE/weight - ! Sample cosine of angle -- fission neutrons are always emitted - ! isotropically. Sometimes in ACE data, fission reactions actually have - ! an angular distribution listed, but for those that do, it's simply just - ! a uniform distribution in mu + ! Sample cosine of angle -- fission neutrons are treated as being emitted + ! isotropically. mu = TWO * prn() - ONE ! Sample azimuthal angle uniformly in [0,2*pi) - phi = TWO*PI*prn() + phi = TWO * PI * prn() bank_array(i) % uvw(1) = mu bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) @@ -256,7 +255,7 @@ contains ! Sample secondary energy distribution for fission reaction and set energy ! in fission bank bank_array(i) % E = & - real(xs % sample_fission_energy(p % g, fission_bank(i) % uvw), 8) + real(xs % sample_fission_energy(p % g, bank_array(i) % uvw), 8) end do ! increment number of bank sites diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index f8fddbc6b3..12c11e2e4e 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -8,36 +8,56 @@ module scattdata_header implicit none + +!=============================================================================== +! JAGGED1D and JAGGED2D is a type which allows for jagged 1-D or 2-D array. +!=============================================================================== + + type :: Jagged2D + real(8), allocatable :: data(:, :) + end type Jagged2D + + type :: Jagged1D + real(8), allocatable :: data(:) + end type Jagged1D + !=============================================================================== ! SCATTDATA contains all the data to describe the scattering energy and ! angular distribution !=============================================================================== type, abstract :: ScattData - ! p0 matrix on its own for sampling energy - real(8), allocatable :: energy(:,:) ! (Gout x Gin) - real(8), allocatable :: mult(:,:) ! (Gout x Gin) - real(8), allocatable :: data(:,:,:) ! (Order/Nmu x Gout x Gin) + ! 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) + type(Jagged1D), allocatable :: mult(:) ! (Gin % data(Gout)) + ! Angular distribution + 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}} contains 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_(this, order, energy, mult, coeffs) + subroutine scattdata_init_(this, mult, coeffs) import ScattData - class(ScattData), intent(inout) :: this ! Object to work on - integer, intent(in) :: order ! Data Order - real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix - 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) import ScattData - class(ScattData), intent(in) :: this ! The ScattData to evaluate + class(ScattData), intent(in) :: this ! Scattering Object to work with integer, intent(in) :: gin ! Incoming Energy Group integer, intent(in) :: gout ! Outgoing Energy Group real(8), intent(in) :: mu ! Angle of interest @@ -47,7 +67,7 @@ module scattdata_header subroutine scattdata_sample_(this, gin, gout, mu, wgt) import ScattData - class(ScattData), intent(in) :: this ! Scattering Object to Use + class(ScattData), intent(in) :: this ! Scattering Object to work with integer, intent(in) :: gin ! Incoming neutron group integer, intent(out) :: gout ! Sampled outgoin group real(8), intent(out) :: mu ! Sampled change in angle @@ -57,30 +77,35 @@ module scattdata_header type, extends(ScattData) :: ScattDataLegendre ! Maximal value for rejection sampling from rectangle - real(8), allocatable :: max_val(:,:) + 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 => scattdatalegendre_init + procedure :: calc_f => scattdatalegendre_calc_f + procedure :: sample => scattdatalegendre_sample end type ScattDataLegendre - type, extends(ScattData) :: ScattDataHistogram - real(8), allocatable :: mu(:) ! Mu bins - real(8) :: dmu ! Mu spacing + type, extends(ScattData) :: ScattDataHistogram + real(8), allocatable :: mu(:) ! Mu bins + real(8) :: dmu ! Mu spacing + ! 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 :: init => scattdatahistogram_init + procedure :: calc_f => scattdatahistogram_calc_f + procedure :: sample => scattdatahistogram_sample + procedure :: get_matrix => scattdatahistogram_get_matrix end type ScattDataHistogram - type, extends(ScattData) :: ScattDataTabular - real(8), allocatable :: mu(:) ! Mu bins - real(8) :: dmu ! Mu spacing - real(8), allocatable :: fmu(:,:,:) ! PDF of f(mu) + type, extends(ScattData) :: ScattDataTabular + real(8), allocatable :: mu(:) ! Mu bins + real(8) :: dmu ! Mu spacing + ! 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 :: init => scattdatatabular_init + procedure :: calc_f => scattdatatabular_calc_f + procedure :: sample => scattdatatabular_sample + procedure :: get_matrix => scattdatatabular_get_matrix end type ScattDataTabular !=============================================================================== @@ -94,189 +119,313 @@ module scattdata_header contains !=============================================================================== -! SCATTDATA_INIT builds the scattdata object +! SCATTDATA*_INIT builds the scattdata object !=============================================================================== 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(in) :: 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 + integer :: groups, gmin, gmax, gin + real(8) :: norm groups = size(energy, dim=1) - allocate(this % energy(groups, groups)) - this % energy = energy - allocate(this % mult(groups, groups)) - this % mult = mult - allocate(this % data(order, groups, groups)) - this % data = ZERO - + allocate(this % gmin(groups)) + allocate(this % gmax(groups)) + 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 + 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)) + this % dist(gin) % data = ZERO + this % gmin(gin) = gmin + this % gmax(gin) = gmax + end do end subroutine scattdata_init - subroutine scattdatalegendre_init(this, order, energy, mult, coeffs) - class(ScattDataLegendre), intent(inout) :: this ! Object to work on - integer, intent(in) :: order ! Data Order - real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix - real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix - real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use + 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 - real(8) :: dmu, mu, f - integer :: imu, Nmu, gout, gin, groups + real(8) :: dmu, mu, f, norm + integer :: imu, Nmu, gout, gin, groups, order + real(8), allocatable :: energy(:, :) + real(8), allocatable :: matrix(:, :, :) - call scattdata_init(this, order, energy, mult) + groups = size(coeffs, dim=3) + order = size(coeffs, dim=1) - this % data = coeffs + ! make a copy of coeffs that we can use to extract data and normalize + allocate(matrix(order, groups, groups)) + matrix (:, :, :)= coeffs - groups = size(this % energy,dim=1) + ! 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) - allocate(this % max_val(groups, groups)) - this % max_val = ZERO - ! Step through the polynomial with fixed number of points to identify - ! the maximal value. - Nmu = 1001 - dmu = TWO / real(Nmu,8) - do imu = 1, Nmu - ! Update mu. Do first and last seperate to avoid float errors - if (imu == 1) then - mu = -ONE - else if (imu == Nmu) then - mu = ONE - end if - mu = -ONE + real(imu - 1,8) * dmu - do gin = 1, groups - do gout = 1, groups - ! Calculate probability - f = this % calc_f(gin,gout,mu) - ! If this is a new max, store it. - if (f > this % max_val(gout,gin)) this % max_val(gout,gin) = f - end do + 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 + if (norm /= ZERO) then + matrix(:, gout, gin) = matrix(:, gout, gin) / norm + end if end do end do - ! Finally, since we may not have caught the exact max, add 10% margin - this % max_val = this % max_val * 1.1_8 + call scattdata_init(this, order, 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) + end do + allocate(this % max_val(gin) % data(this % gmin(gin):this % gmax(gin))) + 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) + do gin = 1, groups + do gout = this % gmin(gin), this % gmax(gin) + do imu = 1, Nmu + ! Update mu. Do first and last seperate to avoid float errors + if (imu == 1) then + mu = -ONE + else if (imu == Nmu) then + mu = ONE + else + mu = -ONE + real(imu - 1, 8) * dmu + end if + ! Calculate probability + f = this % calc_f(gin,gout,mu) + ! If this is a new max, store it. + if (f > this % max_val(gin) % data(gout)) & + this % max_val(gin) % data(gout) = f + end do + ! Finally, since we may not have caught the exact max, add 10% margin + this % max_val(gin) % data(gout) = & + this % max_val(gin) % data(gout) * 1.1_8 + end do + end do end subroutine scattdatalegendre_init - subroutine scattdatahistogram_init(this, order, energy, mult, coeffs) + subroutine scattdatahistogram_init(this, mult, coeffs) class(ScattDataHistogram), intent(inout) :: this ! Object to work on - integer, intent(in) :: order ! Data Order - real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix - 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 + integer :: imu, gin, gout, groups, order real(8) :: norm + real(8), allocatable :: energy(:, :) + real(8), allocatable :: matrix(:, :, :) - groups = size(energy,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 + + ! 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) + + 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 + if (norm /= ZERO) then + matrix(:, gout, gin) = matrix(:, gout, gin) / norm + end if + end do + end do 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 - ! Best to integrate this histogram so we can avoid rejection sampling + ! Integrate this histogram so we can avoid rejection sampling while + ! also saving the original histogram in fmu + allocate(this % fmu(groups)) do gin = 1, groups - do gout = 1, groups - if (energy(gout,gin) > ZERO) then - ! Integrate the histogram - this % data(1,gout,gin) = this % dmu * coeffs(1,gout,gin) - do imu = 2, order - this % data(imu,gout,gin) = this % dmu * coeffs(imu,gout,gin) + & - this % data(imu-1,gout,gin) - end do - ! Now make sure integral norms to zero - norm = this % data(order,gout,gin) - if (norm > ZERO) then - this % data(:,gout,gin) = this % data(:,gout,gin) / norm - end if + allocate(this % fmu(gin) % data(order, & + 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) + ! Integrate the histogram + 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) + end do + + ! Now make sure integral norms to zero + 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 end if end do end do end subroutine scattdatahistogram_init - subroutine scattdatatabular_init(this, order, energy, mult, coeffs) + subroutine scattdatatabular_init(this, mult, coeffs) class(ScattDataTabular), intent(inout) :: this ! Object to work on - integer, intent(in) :: order ! Data Order - real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix - 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 + integer :: imu, gin, gout, groups, order real(8) :: norm - logical :: legendre_flag - integer :: this_order + real(8), allocatable :: energy(:, :) + real(8), allocatable :: matrix(:, :, :) - if (order < 0) then - legendre_flag = .true. - this_order = -1 * order - else - legendre_flag = .false. - this_order = order - end if + groups = size(coeffs, dim=3) + order = size(coeffs, dim=1) - groups = size(energy,dim=1) + ! make a copy of coeffs that we can use to extract data and normalize + allocate(matrix(order, groups, groups)) + matrix(:, :, :) = coeffs - call scattdata_init(this, this_order, energy, mult) - - allocate(this % mu(this_order)) - this % dmu = TWO / real(this_order - 1) - do imu = 1, this_order - 1 - this % mu(imu) = -ONE + real(imu - 1) * this % dmu + ! Build the angular distribution mu values + allocate(this % mu(order)) + 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 end do - this % mu(this_order) = ONE + this % mu(order) = ONE - ! Best to integrate this histogram so we can avoid rejection sampling - allocate(this % fmu(this_order,groups,groups)) + ! Get scattxs + allocate(this % scattxs(groups)) + ! Get this by integrating the scattering distribution over all mu points + ! and then combining over all outgoing groups + ! over all outgoing groups + do gin = 1, groups + norm = ZERO + do gout = 1, groups + do imu = 2, order + 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 + ! Build energy transfer probability matrix from data in matrix do gin = 1, groups do gout = 1, groups - if (energy(gout,gin) > ZERO) then - if (legendre_flag) then - ! Coeffs are legendre coeffs. Need to build f(mu) then integrate - ! and store the integral in this % data - ! Ensure the coeffs are normalized - norm = ONE / coeffs(1,gout,gin) - do imu = 1, this_order - this % fmu(imu,gout,gin) = evaluate_legendre(norm * coeffs(:,gout,gin), this % mu(imu)) - ! Force positivity - if (this % fmu(imu,gout,gin) < ZERO) then - this % fmu(imu,gout,gin) = ZERO - end if - end do - else - ! Coeffs contain f(mu), put in f(mu) to save duplicate. - this % fmu(:,gout,gin) = this % data(:,gout,gin) - end if + norm = ZERO + do imu = 2, order + norm = norm + HALF * this % dmu * & + (matrix(imu - 1, gout, gin) + matrix(imu, gout, gin)) + end do + energy(gout, gin) = norm + end do + end do + call scattdata_init(this, order, energy, mult) - ! Re-normalize fmu for numerical integration issues and in case - ! the negative fix-up introduced un-normalized data - norm = ZERO - do imu = 2, this_order - norm = norm + HALF * this % dmu * (this % fmu(imu-1,gout,gin) + this % fmu(imu,gout,gin)) - end do - if (norm > ZERO) then - this % fmu(:,gout,gin) = this % fmu(:,gout,gin) / norm - end if + ! 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) + ! Coeffs contain f(mu), put in f(mu) as that is where the + ! PDF lives + this % fmu(gin) % data(:, gout) = matrix(:, gout, gin) - ! Now create CDF from fmu with trapezoidal rule - this % data(1,gout,gin) = ZERO - do imu = 2, this_order - 1 - this % data(imu,gout,gin) = this % data(imu-1,gout,gin) + & - HALF * this % dmu * (this % fmu(imu-1,gout,gin) + this % fmu(imu,gout,gin)) - end do - this % data(this_order,gout,gin) = ONE + ! Force positivity + do imu = 1, order + if (this % fmu(gin) % data(imu, gout) < ZERO) then + this % fmu(gin) % data(imu, gout) = ZERO + end if + end do + + ! Re-normalize fmu for numerical integration issues and in case + ! the negative fix-up introduced un-normalized data + norm = ZERO + do imu = 2, order + norm = norm + HALF * this % dmu * & + (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 + end if + + ! Now create CDF from fmu with trapezoidal rule + 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)) + end do + ! Ensure we normalize to 1 still + norm = this % dist(gin) % data(order, gout) + if (norm > ZERO) then + this % dist(gin) % data(:, gout) = & + this % dist(gin) % data(:, gout) / norm end if end do end do - end subroutine scattdatatabular_init !=============================================================================== @@ -285,58 +434,69 @@ contains pure function scattdatalegendre_calc_f(this, gin, gout, mu) result(f) class(ScattDataLegendre), intent(in) :: this ! The ScattData to evaluate - integer, intent(in) :: gin ! Incoming Energy Group - integer, intent(in) :: gout ! Outgoing Energy Group - real(8), intent(in) :: mu ! Angle of interest - real(8) :: f ! Return value of f(mu) - - ! Plug mu in to the legendre expansion and go from there - f = evaluate_legendre(this % data(:, gout, gin), mu) - - end function scattdatalegendre_calc_f - - pure function scattdatahistogram_calc_f(this, gin, gout, mu) result(f) - class(ScattDataHistogram), intent(in) :: this ! The ScattData to evaluate - integer, intent(in) :: gin ! Incoming Energy Group - integer, intent(in) :: gout ! Outgoing Energy Group - real(8), intent(in) :: mu ! Angle of interest - real(8) :: f ! Return value of f(mu) - - integer :: imu - - ! Find mu bin - imu = floor((mu + ONE)/ this % dmu + ONE) - ! Adjust so interpolation works on the last bin if necessary - if (imu == size(this % data, dim=1)) then - imu = imu - 1 - end if - - ! Use histogram interpolation to find f(mu) - f = this % data(imu, gout, gin) - - end function scattdatahistogram_calc_f - - pure function scattdatatabular_calc_f(this, gin, gout, mu) result(f) - class(ScattDataTabular), intent(in) :: this ! The ScattData to evaluate integer, intent(in) :: gin ! Incoming Energy Group integer, intent(in) :: gout ! Outgoing Energy Group real(8), intent(in) :: mu ! Angle of interest real(8) :: f ! Return value of f(mu) + ! Plug mu in to the legendre expansion and go from there + if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then + f = ZERO + else + f = evaluate_legendre(this % dist(gin) % data(:, gout), mu) + end if + + end function scattdatalegendre_calc_f + + pure function scattdatahistogram_calc_f(this, gin, gout, mu) result(f) + class(ScattDataHistogram), intent(in) :: this ! The ScattData to evaluate + integer, intent(in) :: gin ! Incoming Energy Group + integer, intent(in) :: gout ! Outgoing Energy Group + real(8), intent(in) :: mu ! Angle of interest + real(8) :: f ! Return value of f(mu) + + integer :: imu + + if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then + f = ZERO + else + ! Find mu bin + if (mu == ONE) then + imu = size(this % fmu(gin) % data, dim=1) + else + imu = floor((mu + ONE) / this % dmu + ONE) + end if + + f = this % fmu(gin) % data(imu, gout) + end if + + end function scattdatahistogram_calc_f + + pure function scattdatatabular_calc_f(this, gin, gout, mu) result(f) + class(ScattDataTabular), intent(in) :: this ! The ScattData to evaluate + integer, intent(in) :: gin ! Incoming Energy Group + integer, intent(in) :: gout ! Outgoing Energy Group + real(8), intent(in) :: mu ! Angle of interest + real(8) :: f ! Return value of f(mu) + integer :: imu real(8) :: r - ! Find mu bin - imu = floor((mu + ONE)/ this % dmu + ONE) - ! Adjust so interpolation works on the last bin if necessary - if (imu == size(this % data, dim=1)) then - imu = imu - 1 - end if + if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then + f = ZERO + else + ! Find mu bin + if (mu == ONE) then + imu = size(this % fmu(gin) % data, dim=1) - 1 + else + 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 % data(imu, gout, gin) + & - r * this % data(imu + 1, gout, gin) + ! 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) + end if end function scattdatatabular_calc_f @@ -357,24 +517,24 @@ contains integer :: samples xi = prn() - prob = ZERO - gout = 0 + gout = this % gmin(gin) + prob = this % energy(gin) % data(gout) do while (prob < xi) gout = gout + 1 - prob = prob + this % energy(gout,gin) + prob = prob + this % energy(gin) % data(gout) end do - ! Now we can sample mu using the legendre representation of the thisering + ! Now we can sample mu using the legendre representation of the scattering ! kernel in data(1:this % order) - ! Do with rejection sampling + ! Do with rejection sampling from a rectangular bounding box ! Set maximal value - M = this % max_val(gout,gin) + M = this % max_val(gin) % data(gout) 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 @@ -387,7 +547,7 @@ contains end if end do - wgt = wgt * this % mult(gout,gin) + wgt = wgt * this % mult(gin) % data(gout) end subroutine scattdatalegendre_sample @@ -403,26 +563,26 @@ contains integer :: imu xi = prn() - prob = ZERO - gout = 0 + gout = this % gmin(gin) + prob = this % energy(gin) % data(gout) do while (prob < xi) gout = gout + 1 - prob = prob + this % energy(gout,gin) + prob = prob + this % energy(gin) % data(gout) end do xi = prn() - if (xi < this % data(1,gout,gin)) then + if (xi < this % dist(gin) % data(1, gout)) then imu = 1 else - imu = binary_search(this % data(:,gout,gin), & - size(this % data(:,gout,gin)), 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. mu = prn() * this % dmu + this % mu(imu) - wgt = wgt * this % mult(gout,gin) + wgt = wgt * this % mult(gin) % data(gout) end subroutine scattdatahistogram_sample @@ -440,21 +600,21 @@ contains integer :: k, NP xi = prn() - prob = ZERO - gout = 0 + gout = this % gmin(gin) + prob = this % energy(gin) % data(gout) do while (prob < xi) gout = gout + 1 - prob = prob + this % energy(gout,gin) + prob = prob + this % energy(gin) % data(gout) end do ! determine outgoing cosine bin - NP = size(this % data(:,gout,gin)) + NP = size(this % dist(gin) % data(:, gout)) xi = prn() - c_k = this % data(1,gout,gin) + c_k = this % dist(gin) % data(1, gout) do k = 1, NP - 1 - c_k1 = this % data(k+1,gout,gin) + c_k1 = this % dist(gin) % data(k + 1, gout) if (xi < c_k1) exit c_k = c_k1 end do @@ -462,18 +622,19 @@ contains ! check to make sure k is <= NP - 1 k = min(k, NP - 1) - p0 = this % fmu(k,gout,gin) + p0 = this % fmu(gin) % data(k, gout) mu0 = this % mu(k) ! Linear-linear interpolation to find mu value w/in bin. - p1 = this % fmu(k+1,gout,gin) - mu1 = this % mu(k+1) + 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 @@ -482,8 +643,86 @@ contains mu = ONE end if - wgt = wgt * this % mult(gout,gin) + wgt = wgt * this % mult(gin) % data(gout) end subroutine scattdatatabular_sample +!=============================================================================== +! SCATTDATA*_GET_MATRIX Reproduces the original scattering matrix (densely) +! 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 + + integer :: order, groups, gin, gout + + groups = size(this % energy) + ! Set gin and gout for getting the order + order = min(req_order, size(this % dist(1) % data, dim=1)) + + 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 + do gin = 1, groups + do gout = this % gmin(gin), this % gmax(gin) + matrix(:, gout, gin) = this % scattxs(gin) * & + this % energy(gin) % data(gout) * & + this % dist(gin) % data(1:order, gout) + end do + end do + end function 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 + + 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)) + ! 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 + do gout = this % gmin(gin), this % gmax(gin) + matrix(:, gout, gin) = this % scattxs(gin) * & + this % energy(gin) % data(gout) * & + this % fmu(gin) % data(1:order, gout) + end do + end do + end function 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 + + 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)) + ! 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 + do gout = this % gmin(gin), this % gmax(gin) + matrix(:, gout, gin) = 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 module scattdata_header diff --git a/src/state_point.F90 b/src/state_point.F90 index 10979ee9e0..3970a7331c 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -296,7 +296,11 @@ contains NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins if (tally % nuclide_bins(j) > 0) then ! Get index in cross section listings for this nuclide - i_list = nuclides(tally % nuclide_bins(j)) % listing + if (run_CE) then + i_list = nuclides(tally % nuclide_bins(j)) % listing + else + i_list = nuclides_MG(tally % nuclide_bins(j)) % obj % listing + end if ! Determine position of . in alias string (e.g. "U-235.71c"). If ! no . is found, just use the entire string. diff --git a/src/summary.F90 b/src/summary.F90 index 35c62ae039..4502058cad 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -653,7 +653,11 @@ contains allocate(str_array(t%n_nuclide_bins)) NUCLIDE_LOOP: do j = 1, t%n_nuclide_bins if (t%nuclide_bins(j) > 0) then - i_list = nuclides(t%nuclide_bins(j))%listing + if (run_CE) then + i_list = nuclides(t % nuclide_bins(j)) % listing + else + i_list = nuclides_MG(t % nuclide_bins(j)) % obj % listing + end if i_xs = index(xs_listings(i_list)%alias, '.') if (i_xs > 0) then str_array(j) = xs_listings(i_list)%alias(1:i_xs - 1) diff --git a/src/tally.F90 b/src/tally.F90 index 1f0d49772b..c4eaf30c8c 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -26,8 +26,9 @@ module tally !$omp threadprivate(position) - procedure(score_general_), pointer :: score_general => null() - procedure(get_scoring_bins_), pointer :: get_scoring_bins => null() + procedure(score_general_), pointer :: score_general => null() + procedure(score_analog_tally_), pointer :: score_analog_tally => null() + procedure(get_scoring_bins_), pointer :: get_scoring_bins => null() abstract interface subroutine score_general_(p, t, start_index, filter_index, i_nuclide, & @@ -43,6 +44,11 @@ module tally real(8), intent(in) :: atom_density ! atom/b-cm end subroutine score_general_ + subroutine score_analog_tally_(p) + import Particle + type(Particle), intent(in) :: p + end subroutine score_analog_tally_ + subroutine get_scoring_bins_(p, i_tally, found_bin) import Particle type(Particle), intent(in) :: p @@ -60,11 +66,13 @@ contains subroutine init_tally_routines() if (run_CE) then - score_general => score_general_ce - get_scoring_bins => get_scoring_bins_ce + score_general => score_general_ce + score_analog_tally => score_analog_tally_ce + get_scoring_bins => get_scoring_bins_ce else - score_general => score_general_mg - get_scoring_bins => get_scoring_bins_mg + score_general => score_general_mg + score_analog_tally => score_analog_tally_mg + get_scoring_bins => get_scoring_bins_mg end if end subroutine init_tally_routines @@ -801,16 +809,48 @@ contains integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: score_index ! scoring bin index real(8) :: score ! analog tally score - real(8) :: macro_total ! material macro total xs - real(8) :: macro_scatt ! material macro scatt xs - real(8) :: micro_abs ! nuclidic microscopic abs real(8) :: p_uvw(3) ! Particle's current uvw + integer :: p_g ! Particle group to use for getting info + ! to tally with. + class(Mgxs), pointer :: matxs + class(Mgxs), pointer :: nucxs - ! Set the direction, if needed for nuclidic data, so that nuc % get_xs - ! knows wihch direction it should be using for direction-dependent - ! mgxs - if (i_nuclide > 0) then + ! Set the direction and group to use with get_xs + ! this only depends on if we + if (t % estimator == ESTIMATOR_ANALOG .or. & + t % estimator == ESTIMATOR_COLLISION) then + if (survival_biasing) then + ! Then we either are alive and had a scatter (and so g changed), + ! or are dead and g did not change + if (p % alive) then + p_uvw = p % last_uvw + p_g = p % last_g + else + p_uvw = p % coord(p % n_coord) % uvw + p_g = p % g + end if + else if (p % event == EVENT_SCATTER) then + ! Then the energy group has been changed by the scattering routine + ! meaning gin is now in p % last_g + p_uvw = p % last_uvw + p_g = p % last_g + else + ! No scatter, no change in g. + p_uvw = p % coord(p % n_coord) % uvw + p_g = p % g + end if + else + ! No actual collision so g has not changed. p_uvw = p % coord(p % n_coord) % uvw + p_g = p % g + end if + + ! To significantly reduce de-referencing, point matxs to the + ! macroscopic Mgxs for the material of interest + matxs => macro_xs(p % material) % obj + ! Do same for nucxs, point it to the microscopic nuclide data of interest + if (i_nuclide > 0) then + nucxs => nuclides_MG(i_nuclide) % obj end if i = 0 @@ -861,13 +901,16 @@ contains else score = p % last_wgt end if + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('total', p_g, UVW=p_uvw) / & + matxs % get_xs('total', p_g, UVW=p_uvw) + end if else if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = nuc % get_xs(p % g, 'total', UVW=p_uvw) * & - atom_density * flux - end associate + score = nucxs % get_xs('total', p_g, UVW=p_uvw) * & + atom_density * flux else score = material_xs % total * flux end if @@ -875,7 +918,8 @@ contains case (SCORE_INVERSE_VELOCITY) - if (t % estimator == ESTIMATOR_ANALOG) then + if (t % estimator == ESTIMATOR_ANALOG .or. & + t % estimator == ESTIMATOR_COLLISION) then ! All events score to an inverse velocity bin. We actually use a ! collision estimator in place of an analog one since there is no way ! to count 'events' exactly for the inverse velocity @@ -886,150 +930,94 @@ contains else score = p % last_wgt end if - score = score * inverse_velocities(p % last_g) + score = score * inverse_velocities(p_g) / material_xs % total else ! For inverse velocity, we need no cross section - score = score * inverse_velocities(p % g) + score = flux * inverse_velocities(p_g) end if - case (SCORE_SCATTER, SCORE_SCATTER_N) + case (SCORE_SCATTER, SCORE_SCATTER_N, SCORE_SCATTER_PN, SCORE_SCATTER_YN) if (t % estimator == ESTIMATOR_ANALOG) then ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + if (p % event /= EVENT_SCATTER) then + if (score_bin == SCORE_SCATTER_PN) then + i = i + t % moment_order(i) + else if (score_bin == SCORE_SCATTER_YN) then + i = i + (t % moment_order(i) + 1)**2 - 1 + end if + cycle SCORE_LOOP + end if + ! Since only scattering events make it here, again we can use ! the weight entering the collision as the estimator for the ! reaction rate score = p % last_wgt - else - ! Note SCORE_SCATTER_N not available for tracklength/collision. + ! Since we transport based on material data, the angle selected + ! was not selected from the f(mu) for the nuclide. Therefore + ! adjust the score by the actual probability for that nuclide. if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = nuc % get_xs(p % g, 'scatter', UVW=p_uvw) * & - atom_density * flux - end associate + score = score * atom_density * & + nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, & + UVW=p_uvw, MU=p % mu) / & + matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, & + UVW=p_uvw, MU=p % mu) + end if + + else + ! Note SCORE_SCATTER_*N not available for tracklength/collision. + if (i_nuclide > 0) then + score = atom_density * flux * & + nucxs % get_xs('scatter/mult', p_g, UVW=p_uvw) else - ! Get the scattering x/s (stored in % elastic) - score = material_xs % elastic * flux + ! Get the scattering x/s and take away + ! the multiplication baked in to sigS + score = flux * & + matxs % get_xs('scatter/mult', p_g, UVW=p_uvw) end if end if - if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, & - p % last_uvw, p % mu) - end associate + + case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N, SCORE_NU_SCATTER_PN, & + SCORE_NU_SCATTER_YN) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + if (score_bin == SCORE_NU_SCATTER_PN) then + i = i + t % moment_order(i) + else if (score_bin == SCORE_NU_SCATTER_YN) then + i = i + (t % moment_order(i) + 1)**2 - 1 + end if + cycle SCORE_LOOP + end if + + ! For scattering production, we need to use the pre-collision + ! weight times the multiplicity as the estimate for the number of + ! neutrons exiting a reaction with neutrons in the exit channel + score = p % wgt + + ! Since we transport based on material data, the angle selected + ! was not selected from the f(mu) for the nuclide. Therefore + ! adjust the score by the actual probability for that nuclide. + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('scatter*f_mu', p % last_g, p % g, & + UVW=p_uvw, MU=p % mu) / & + matxs % get_xs('scatter*f_mu', p % last_g, p % g, & + UVW=p_uvw, MU=p % mu) + end if + else - score = score / & - macro_xs(p % material) % obj % get_xs(p % g, 'mult', & - p % last_g, & - p % last_uvw) - end if - - - case (SCORE_SCATTER_PN) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - i = i + t % moment_order(i) - cycle SCORE_LOOP - end if - ! Since only scattering events make it here, again we can use - ! the weight entering the collision as the estimator for the - ! reaction rate - score = p % last_wgt - - if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, & - p % last_uvw, p % mu) - end associate - else - score = score / & - macro_xs(p % material) % obj % get_xs(p % g, 'mult', & - p % last_g, & - p % last_uvw) - end if - - - case (SCORE_SCATTER_YN) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - i = i + (t % moment_order(i) + 1)**2 - 1 - cycle SCORE_LOOP - end if - ! Since only scattering events make it here, again we can use - ! the weight entering the collision as the estimator for the - ! reaction rate - score = p % last_wgt - - if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, & - p % last_uvw, p % mu) - end associate - else - score = score / & - macro_xs(p % material) % obj % get_xs(p % g, 'mult', & - p % last_g, & - p % last_uvw) - end if - - - case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! For scattering production, we need to use the pre-collision - ! weight times the multiplicity as the estimate for the number of - ! neutrons exiting a reaction with neutrons in the exit channel - score = p % wgt - if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, & - p % last_uvw, p % mu) - end associate - end if - - - case (SCORE_NU_SCATTER_PN) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - i = i + t % moment_order(i) - cycle SCORE_LOOP - end if - ! For scattering production, we need to use the pre-collision - ! weight times the multiplicity as the estimate for the number of - ! neutrons exiting a reaction with neutrons in the exit channel - score = p % wgt - if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, & - p % last_uvw, p % mu) - end associate - end if - - - case (SCORE_NU_SCATTER_YN) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - i = i + (t % moment_order(i) + 1)**2 - 1 - cycle SCORE_LOOP - end if - ! For scattering production, we need to use the pre-collision - ! weight times the multiplicity as the estimate for the number of - ! neutrons exiting a reaction with neutrons in the exit channel - score = p % wgt - if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, & - p % last_uvw, p % mu) - end associate + ! Note SCORE_NU_SCATTER_*N not available for tracklength/collision. + if (i_nuclide > 0) then + score = nucxs % get_xs('scatter', p_g, UVW=p_uvw) * & + atom_density * flux + else + ! Get the scattering x/s, which includes multiplication + score = matxs % get_xs('scatter', p_g, UVW=p_uvw) * flux + end if end if @@ -1037,24 +1025,14 @@ contains ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! get material macros - macro_total = material_xs % total - macro_scatt = material_xs % elastic ! Score total rate - p1 scatter rate Note estimator needs to be ! adjusted since tallying is only occuring when a scatter has ! happened. Effectively this means multiplying the estimator by ! total/scatter macro - score = (macro_total - p % mu * macro_scatt) * (ONE / macro_scatt) - - - case (SCORE_N_1N) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! Skip any events where weight of particle changed - if (p % wgt /= p % last_wgt) cycle SCORE_LOOP - ! All events that reach this point are (n,1n) reactions - score = p % last_wgt + score = (material_xs % total - p % mu * material_xs % elastic) + if (material_xs % elastic /= ZERO) then + score = score / material_xs % elastic + end if case (SCORE_ABSORPTION) @@ -1070,13 +1048,15 @@ contains ! can just use the particle's weight entering the collision score = p % last_wgt end if - + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % 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 - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = nuc % get_xs(p % g, 'absorption', UVW=p_uvw) & - * atom_density * flux - end associate + score = nucxs % get_xs('absorption', p_g, UVW=p_uvw) * & + atom_density * flux else score = material_xs % absorption * flux end if @@ -1089,38 +1069,30 @@ contains ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in ! fission - associate (nuc => nuclides_MG(i_nuclide) % obj) - micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p_uvw) - if (micro_abs > ZERO) then - score = p % absorb_wgt * & - nuc % get_xs(p % g, 'fission', UVW=p_uvw) & - / micro_abs - else - score = ZERO - end if - end associate + score = p % absorb_wgt else ! Skip any non-absorption events if (p % event == EVENT_SCATTER) cycle SCORE_LOOP ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for the ! fission reaction rate - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = p % last_wgt & - * nuc % get_xs(p % g, 'fission', UVW=p_uvw) & - / nuc % get_xs(p % g, 'absorption', UVW=p_uvw) - end associate + score = p % last_wgt end if - + if (i_nuclide > 0) then + score = score * atom_density * & + 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('absorption', p_g, UVW=p_uvw) + end if else if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = nuc % get_xs(p % g, 'fission', UVW=p_uvw) * & - atom_density * flux - end associate + score = nucxs % get_xs('fission', p_g, UVW=p_uvw) * & + atom_density * flux else - score = flux * macro_xs(p % material) % obj % get_xs(p % g, & - 'fission', UVW=p_uvw) + score = flux * material_xs % fission end if end if @@ -1135,7 +1107,8 @@ contains ! neutrons were emitted with different energies, multiple ! outgoing energy bins may have been scored to. The following ! logic treats this special case and results to multiple bins - call score_fission_eout_mg(p, t, score_index) + call score_fission_eout_mg(p, t, score_index, i_nuclide, & + atom_density) cycle SCORE_LOOP end if end if @@ -1143,16 +1116,16 @@ contains ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in ! nu-fission - associate (nuc => nuclides_MG(i_nuclide) % obj) - micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p_uvw) - if (micro_abs > ZERO) then - score = p % absorb_wgt * & - nuc % get_xs(p % g, 'fission', UVW=p_uvw) / & - micro_abs - else - score = ZERO - end if - end associate + score = p % absorb_wgt + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) + else + score = score * & + matxs % get_xs('nu_fission', p_g, UVW=p_uvw) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) + end if else ! Skip any non-fission events if (.not. p % fission) cycle SCORE_LOOP @@ -1162,14 +1135,17 @@ contains ! bank. Since this was weighted by 1/keff, we multiply by keff ! to get the proper score. score = keff * p % wgt_bank + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('fission', p_g, UVW=p_uvw) / & + matxs % get_xs('fission', p_g, UVW=p_uvw) + end if end if else if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = nuc % get_xs(p % g, 'nu_fission', UVW=p_uvw) & - * atom_density * flux - end associate + score = nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) * & + atom_density * flux else score = material_xs % nu_fission * flux end if @@ -1177,43 +1153,36 @@ contains case (SCORE_KAPPA_FISSION) - ! Determine kappa-fission cross section - score = ZERO if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in - ! fission scale by kappa-fission - associate (nuc => nuclides_MG(i_nuclide) % obj) - micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p_uvw) - if (micro_abs > ZERO) then - score = p % absorb_wgt * & - nuc % get_xs(p % g, 'k_fission', UVW=p_uvw) / & - micro_abs - end if - end associate + ! fission + score = p % absorb_wgt else ! Skip any non-absorption events if (p % event == EVENT_SCATTER) cycle SCORE_LOOP ! All fission events will contribute, so again we can use - ! particle's weight entering the collision as the estimate for - ! the fission energy production rate - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = p % last_wgt * & - nuc % get_xs(p % g, 'k_fission', UVW=p_uvw) / & - nuc % get_xs(p % g, 'absorption', UVW=p_uvw) - end associate + ! particle's weight entering the collision as the estimate for the + ! fission reaction rate + score = p % last_wgt + end if + 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) + else + score = score * & + matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) end if - else if (i_nuclide > 0) then - associate (nuc => nuclides_MG(i_nuclide) % obj) - score = nuc % get_xs(p % g, 'k_fission', UVW=p_uvw) & - * atom_density * flux - end associate + score = flux * atom_density * & + nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) else - score = flux * macro_xs(p % material) % obj % get_xs(p % g, & - 'k_fission', UVW=p_uvw) + score = flux * matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) + end if end if @@ -1230,6 +1199,8 @@ contains score, i) end do SCORE_LOOP + + nullify(matxs,nucxs) end subroutine score_general_mg !=============================================================================== @@ -1404,7 +1375,7 @@ contains ! triggered at every collision, not every event !=============================================================================== - subroutine score_analog_tally(p) + subroutine score_analog_tally_ce(p) type(Particle), intent(in) :: p @@ -1414,17 +1385,9 @@ contains ! position during the loop integer :: filter_index ! single index for single bin integer :: i_nuclide ! index in nuclides array - real(8) :: last_wgt ! pre-collision particle weight - real(8) :: wgt ! post-collision particle weight - real(8) :: mu ! cosine of angle of collision logical :: found_bin ! scoring bin found? type(TallyObject), pointer :: t - ! Copy particle's pre- and post-collision weight and angle - last_wgt = p % last_wgt - wgt = p % wgt - mu = p % mu - ! A loop over all tallies is necessary because we need to simultaneously ! determine different filter bins for the same tally in order to score to it @@ -1503,7 +1466,87 @@ contains ! Reset tally map positioning position = 0 - end subroutine score_analog_tally + end subroutine score_analog_tally_ce + + subroutine score_analog_tally_mg(p) + + type(Particle), intent(in) :: p + + integer :: i, m + integer :: i_tally + integer :: k ! loop index for nuclide bins + ! position during the loop + integer :: filter_index ! single index for single bin + integer :: i_nuclide ! index in nuclides array + logical :: found_bin ! scoring bin found? + type(TallyObject), pointer :: t + type(Material), pointer :: mat + real(8) :: atom_density + + ! A loop over all tallies is necessary because we need to simultaneously + ! determine different filter bins for the same tally in order to score to it + + TALLY_LOOP: do i = 1, active_analog_tallies % size() + ! Get index of tally and pointer to tally + i_tally = active_analog_tallies % get_item(i) + t => tallies(i_tally) + + ! Get pointer to current material. We need this in order to determine what + ! nuclides are in the material + mat => materials(p % material) + + ! ======================================================================= + ! DETERMINE SCORING BIN COMBINATION + + call get_scoring_bins(p, i_tally, found_bin) + if (.not. found_bin) cycle + + ! ======================================================================= + ! CALCULATE RESULTS AND ACCUMULATE TALLY + + ! If we have made it here, we have a scoring combination of bins for this + ! tally -- now we need to determine where in the results array we should + ! be accumulating the tally values + + ! Determine scoring index for this filter combination + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 + + ! Check for nuclide bins + k = 0 + NUCLIDE_LOOP: do while (k < t % n_nuclide_bins) + + ! Increment the index in the list of nuclide bins + k = k + 1 + + i_nuclide = t % nuclide_bins(k) + + ! Check to see if this nuclide was in the material of our collision. + do m = 1, mat % n_nuclides + if (mat % nuclide(m) == i_nuclide) then + atom_density = mat % atom_density(m) + exit + end if + end do + + ! Determine score for each bin + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, ZERO) + + end do NUCLIDE_LOOP + + ! If the user has specified that we can assume all tallies are spatially + ! separate, this implies that once a tally has been scored to, we needn't + ! check the others. This cuts down on overhead when there are many + ! tallies specified + + if (assume_separate) exit TALLY_LOOP + + end do TALLY_LOOP + + ! Reset tally map positioning + position = 0 + + end subroutine score_analog_tally_mg !=============================================================================== ! SCORE_FISSION_EOUT handles a special case where we need to store neutron @@ -1566,10 +1609,12 @@ contains end subroutine score_fission_eout_ce - subroutine score_fission_eout_mg(p, t, i_score) + subroutine score_fission_eout_mg(p, t, i_score, i_nuclide, atom_density) type(Particle), intent(in) :: p type(TallyObject), intent(inout) :: t - integer, intent(in) :: i_score ! index for score + integer, intent(in) :: i_score ! index for score + integer, intent(in) :: i_nuclide ! index for nuclide + real(8), intent(in) :: atom_density integer :: i ! index of outgoing energy filter integer :: n ! number of energies on filter @@ -1578,6 +1623,7 @@ contains integer :: i_filter ! index for matching filter bin combination real(8) :: score ! actual score integer :: gout ! energy group of fission bank site + integer :: gin ! energy group of incident particle real(8) :: E_out ! save original outgoing energy bin and score index @@ -1596,6 +1642,18 @@ contains do k = 1, p % n_bank ! determine score based on bank site weight and keff score = keff * fission_bank(n_bank - p % n_bank + k) % wgt + if (i_nuclide > 0) then + if (survival_biasing) then + gin = p % g + else + gin = p % last_g + end if + score = score * atom_density * & + nuclides_MG(i_nuclide) % obj % get_xs('fission', gin, & + UVW=p % last_uvw) / & + macro_xs(p % material) % obj % get_xs('fission', gin, & + UVW=p % last_uvw) + end if if (t % energyout_matches_groups) then ! determine outgoing energy from fission bank @@ -1605,7 +1663,7 @@ contains matching_bins(i) = gout else ! determine outgoing energy from fission bank - E_out = fission_bank(n_bank - p % n_bank + k) % E + E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) % E)) ! check if outgoing energy is within specified range on filter if (E_out < t % filters(i) % real_bins(1) .or. & @@ -2601,7 +2659,6 @@ contains end if end if - case (FILTER_ENERGYOUT) if (t % energyout_matches_groups) then ! Since all groups are filters, the filter bin is the group @@ -2625,7 +2682,6 @@ contains end if end if - case (FILTER_MU) ! determine mu bin n = t % filters(i) % n_bins diff --git a/src/tracking.F90 b/src/tracking.F90 index e634112bc6..d52b09a75d 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -7,7 +7,6 @@ module tracking cross_lattice, check_cell_overlap use geometry_header, only: Universe, BASE_UNIVERSE use global - use macroxs_header, only: MacroXS use output, only: write_message use particle_header, only: LocalCoord, Particle use physics, only: collision @@ -98,6 +97,7 @@ contains material_xs % total = ZERO material_xs % elastic = ZERO material_xs % absorption = ZERO + material_xs % fission = ZERO material_xs % nu_fission = ZERO end if end if diff --git a/tests/test_mg_basic/results_true.dat b/tests/test_mg_basic/results_true.dat index 35f3e73d40..55c2af8136 100644 --- a/tests/test_mg_basic/results_true.dat +++ b/tests/test_mg_basic/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.045320E+00 5.851680E-02 +1.033731E+00 4.974463E-02 diff --git a/tests/test_mg_max_order/inputs_true.dat b/tests/test_mg_max_order/inputs_true.dat index 913f8200fb..63bdaab03e 100644 --- a/tests/test_mg_max_order/inputs_true.dat +++ b/tests/test_mg_max_order/inputs_true.dat @@ -1 +1 @@ -7f7465abaf559b3ef56cb6b0f28050c12f392f55db33dc5d2cefc14b92beb2c9068834c05273e51323d3516643e8a385e4c177a7a471678c961808d19055a30f \ No newline at end of file +60a35864ad71646309d7f1687ba0826d4d53a5b2e8babf73614362645205484bad3c0e7bf605ec0b11cadf58474b2e3d0a97bf2d9297f9118682c37ff0269afd \ No newline at end of file diff --git a/tests/test_mg_max_order/results_true.dat b/tests/test_mg_max_order/results_true.dat index 1b23005602..f75c1300ab 100644 --- a/tests/test_mg_max_order/results_true.dat +++ b/tests/test_mg_max_order/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.083030E+00 1.855038E-02 +1.055274E+00 1.715904E-02 diff --git a/tests/test_mg_max_order/test_mg_max_order.py b/tests/test_mg_max_order/test_mg_max_order.py index 7f59572cff..088f6914ba 100644 --- a/tests/test_mg_max_order/test_mg_max_order.py +++ b/tests/test_mg_max_order/test_mg_max_order.py @@ -76,9 +76,10 @@ class MGMaxOrderTestHarness(PyAPITestHarness): self._input_set = MGNuclideInputSet() def _build_inputs(self): - super(MGMaxOrderTestHarness, self)._build_inputs() # Set P1 scattering self._input_set.settings.max_order = 1 + # Call standard input build + super(MGMaxOrderTestHarness, self)._build_inputs() if __name__ == '__main__': harness = MGMaxOrderTestHarness('statepoint.10.*', False, mg=True) diff --git a/tests/test_mg_nuclide/inputs_true.dat b/tests/test_mg_nuclide/inputs_true.dat index 32a7773c1e..e0af3352b2 100644 --- a/tests/test_mg_nuclide/inputs_true.dat +++ b/tests/test_mg_nuclide/inputs_true.dat @@ -1 +1 @@ -825dee3ca35d48788f1a4d5364789bbd83b36e33af9a990da758dd73c3bfcbee14bce2a41e6c80e0147f45575e59078653c8dfa8590cd361c09f19c26dc8c88e \ No newline at end of file +0efba3dd7882fdd38756d0a8f01ff00d7a1abdaab6430b3f090f3339e552448453bbb733852b6bd6ff09608d923c282f168320f942fc2eb3a45610873c588734 \ No newline at end of file diff --git a/tests/test_mg_nuclide/results_true.dat b/tests/test_mg_nuclide/results_true.dat index 5b60cef222..55c2af8136 100644 --- a/tests/test_mg_nuclide/results_true.dat +++ b/tests/test_mg_nuclide/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.380785E-01 5.556526E-03 +1.033731E+00 4.974463E-02 diff --git a/tests/test_mg_nuclide/test_mg_nuclide.py b/tests/test_mg_nuclide/test_mg_nuclide.py index 866840ddf1..89694eea25 100644 --- a/tests/test_mg_nuclide/test_mg_nuclide.py +++ b/tests/test_mg_nuclide/test_mg_nuclide.py @@ -12,16 +12,16 @@ class MGNuclideInputSet(MGInputSet): # Define materials needed for 1D/1G slab problem # This time do using nuclide, not macroscopic uo2 = openmc.Material(name='UO2', material_id=1) - uo2.set_density('g/cm3', 1.0) + uo2.set_density('sum', 1.0) uo2.add_nuclide("uo2_iso", 1.0) clad = openmc.Material(name='Clad', material_id=2) - clad.set_density('g/cm3', 1.0) + clad.set_density('sum', 1.0) clad.add_nuclide("clad_ang_mu", 1.0) - water_data = openmc.Nuclide('lwtr_iso_mu', '71c') + # water_data = openmc.Nuclide('lwtr_iso_mu', '71c') water = openmc.Material(name='LWTR', material_id=3) - water.set_density('g/cm3', 1.0) + water.set_density('sum', 1.0) water.add_nuclide("lwtr_iso_mu", 1.0) # Define the materials file. diff --git a/tests/test_mg_tallies/results_true.dat b/tests/test_mg_tallies/results_true.dat index 0cb47a712d..debbfa537c 100644 --- a/tests/test_mg_tallies/results_true.dat +++ b/tests/test_mg_tallies/results_true.dat @@ -1,86 +1,86 @@ k-combined: -1.045320E+00 5.851680E-02 +1.033731E+00 4.974463E-02 tally 1: -2.286064E+00 -1.057353E+00 -6.503987E-02 -8.851627E-04 -3.376363E+00 -2.323627E+00 -2.733240E-02 -1.607534E-04 -6.776283E-02 -9.880704E-04 -2.391658E+00 -1.201477E+00 -7.241106E-02 -1.103780E-03 -3.614949E+00 -2.730438E+00 -3.146867E-02 -2.110753E-04 -7.801752E-02 -1.297373E-03 -2.762725E+00 -1.705088E+00 -8.684520E-02 -1.654173E-03 -4.172232E+00 -3.847245E+00 -3.834947E-02 -3.212662E-04 -9.507651E-02 -1.974662E-03 -2.802290E+00 -1.773339E+00 -8.347451E-02 -1.574952E-03 -4.206829E+00 -3.971225E+00 -3.593646E-02 -2.967708E-04 -8.909414E-02 -1.824101E-03 -2.383708E+00 -1.176784E+00 -7.337273E-02 -1.097240E-03 -3.624903E+00 -2.690697E+00 -3.213890E-02 -2.139948E-04 -7.967917E-02 -1.315319E-03 -2.398216E+00 -1.234567E+00 -6.905889E-02 -9.879327E-04 -3.479138E+00 -2.538252E+00 -2.911091E-02 -1.750648E-04 -7.217215E-02 -1.076035E-03 -2.563998E+00 -1.354089E+00 -7.357381E-02 -1.097086E-03 -3.753156E+00 -2.867475E+00 -3.097034E-02 -1.948794E-04 -7.678206E-02 -1.197826E-03 -2.293243E+00 -1.172767E+00 -6.702582E-02 -9.267762E-04 -3.407144E+00 -2.469472E+00 -2.857514E-02 -1.688581E-04 -7.084385E-02 -1.037886E-03 +3.163666E+00 +2.165097E+00 +9.964133E-02 +2.052446E-03 +4.861844E+00 +4.978206E+00 +4.417216E-02 +4.040216E-04 +1.095122E-01 +2.483317E-03 +3.324437E+00 +2.299850E+00 +9.574329E-02 +1.968438E-03 +4.881821E+00 +4.943976E+00 +4.041706E-02 +3.650255E-04 +1.002025E-01 +2.243628E-03 +3.199995E+00 +2.091126E+00 +8.859707E-02 +1.592091E-03 +4.671522E+00 +4.439001E+00 +3.660515E-02 +2.728602E-04 +9.075197E-02 +1.677135E-03 +2.910284E+00 +1.723356E+00 +9.207508E-02 +1.744614E-03 +4.421737E+00 +3.979362E+00 +4.080063E-02 +3.481887E-04 +1.011535E-01 +2.140141E-03 +2.506574E+00 +1.326705E+00 +8.637880E-02 +1.598941E-03 +3.920683E+00 +3.263607E+00 +3.978214E-02 +3.424111E-04 +9.862841E-02 +2.104629E-03 +2.951103E+00 +1.826551E+00 +8.748324E-02 +1.648479E-03 +4.309848E+00 +3.903067E+00 +3.741466E-02 +3.141687E-04 +9.275891E-02 +1.931037E-03 +3.048521E+00 +2.007251E+00 +9.483162E-02 +1.935789E-03 +4.599534E+00 +4.527383E+00 +4.168244E-02 +3.813062E-04 +1.033396E-01 +2.343697E-03 +2.982958E+00 +1.966657E+00 +9.896454E-02 +2.016812E-03 +4.645921E+00 +4.599347E+00 +4.489740E-02 +4.108900E-04 +1.113102E-01 +2.525534E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -171,86 +171,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.604127E+00 -1.442914E+00 -7.142299E-02 -1.083324E-03 -3.786547E+00 -2.990260E+00 -2.935394E-02 -1.905369E-04 -7.277467E-02 -1.171135E-03 -2.457755E+00 -1.228862E+00 -6.655630E-02 -9.411086E-04 -3.528688E+00 -2.561047E+00 -2.709124E-02 -1.640319E-04 -6.716494E-02 -1.008221E-03 -2.450846E+00 -1.295337E+00 -6.779278E-02 -9.992248E-04 -3.519409E+00 -2.693620E+00 -2.793186E-02 -1.714064E-04 -6.924902E-02 -1.053549E-03 -2.469234E+00 -1.300419E+00 -7.347034E-02 -1.175743E-03 -3.675903E+00 -2.880386E+00 -3.155996E-02 -2.200036E-04 -7.824386E-02 -1.352252E-03 -2.576106E+00 -1.365945E+00 -7.241428E-02 -1.090052E-03 -3.719498E+00 -2.861357E+00 -3.008961E-02 -1.906170E-04 -7.459854E-02 -1.171627E-03 -2.503651E+00 -1.290812E+00 -7.432507E-02 -1.132918E-03 -3.702398E+00 -2.813425E+00 -3.186491E-02 -2.091477E-04 -7.899991E-02 -1.285525E-03 -2.395349E+00 -1.202098E+00 -7.095925E-02 -1.051988E-03 -3.559388E+00 -2.628753E+00 -3.041477E-02 -1.959992E-04 -7.540470E-02 -1.204709E-03 -1.910174E+00 -7.331782E-01 -6.366813E-02 -8.166135E-04 -2.966770E+00 -1.762376E+00 -2.893278E-02 -1.712979E-04 -7.173053E-02 -1.052882E-03 +2.806322E+00 +1.650116E+00 +8.677689E-02 +1.577748E-03 +4.271173E+00 +3.787053E+00 +3.812147E-02 +3.107861E-04 +9.451124E-02 +1.910246E-03 +2.767101E+00 +1.574082E+00 +8.523508E-02 +1.523190E-03 +4.143067E+00 +3.546386E+00 +3.725690E-02 +2.974028E-04 +9.236779E-02 +1.827986E-03 +2.767691E+00 +1.557585E+00 +7.816996E-02 +1.254835E-03 +4.015666E+00 +3.289915E+00 +3.260527E-02 +2.211412E-04 +8.083542E-02 +1.359243E-03 +2.734236E+00 +1.510579E+00 +8.505695E-02 +1.474976E-03 +4.114798E+00 +3.424404E+00 +3.737976E-02 +2.870194E-04 +9.267239E-02 +1.764164E-03 +2.421441E+00 +1.195490E+00 +8.292263E-02 +1.384777E-03 +3.817448E+00 +2.942804E+00 +3.814207E-02 +2.931372E-04 +9.456230E-02 +1.801767E-03 +2.724650E+00 +1.519118E+00 +9.238437E-02 +1.725646E-03 +4.305602E+00 +3.753830E+00 +4.235399E-02 +3.657508E-04 +1.050046E-01 +2.248086E-03 +2.703678E+00 +1.593893E+00 +8.744732E-02 +1.597761E-03 +4.110489E+00 +3.592877E+00 +3.913038E-02 +3.184391E-04 +9.701255E-02 +1.957285E-03 +2.707705E+00 +1.792230E+00 +8.754036E-02 +1.760104E-03 +4.203619E+00 +4.162497E+00 +3.927598E-02 +3.456207E-04 +9.737351E-02 +2.124357E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -341,86 +341,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.593479E+00 -1.421618E+00 -7.731733E-02 -1.249860E-03 -3.880469E+00 -3.167768E+00 -3.330331E-02 -2.319171E-04 -8.256599E-02 -1.425478E-03 -2.563470E+00 -1.449702E+00 -7.411043E-02 -1.158986E-03 -3.692974E+00 -2.935486E+00 -3.125554E-02 -2.061246E-04 -7.748913E-02 -1.266944E-03 -2.657580E+00 -1.493736E+00 -7.749682E-02 -1.277831E-03 -3.892576E+00 -3.190262E+00 -3.289903E-02 -2.334764E-04 -8.156370E-02 -1.435062E-03 -2.701357E+00 -1.517431E+00 -8.804410E-02 -1.639744E-03 -4.161224E+00 -3.606273E+00 -3.959674E-02 -3.345554E-04 -9.816875E-02 -2.056343E-03 -2.745200E+00 -1.523378E+00 -7.925308E-02 -1.270874E-03 -3.962649E+00 -3.165021E+00 -3.340000E-02 -2.283904E-04 -8.280570E-02 -1.403801E-03 -2.678775E+00 -1.492394E+00 -8.646776E-02 -1.554919E-03 -4.147456E+00 -3.551060E+00 -3.875994E-02 -3.125295E-04 -9.609415E-02 -1.920962E-03 -2.678888E+00 -1.497390E+00 -7.616647E-02 -1.203778E-03 -3.868905E+00 -3.108718E+00 -3.186093E-02 -2.136540E-04 -7.899003E-02 -1.313223E-03 -2.189117E+00 -9.756340E-01 -6.824988E-02 -9.549688E-04 -3.259136E+00 -2.150849E+00 -2.997262E-02 -1.883650E-04 -7.430850E-02 -1.157785E-03 +2.333969E+00 +1.124588E+00 +6.982134E-02 +1.014762E-03 +3.448351E+00 +2.435965E+00 +3.004987E-02 +1.926319E-04 +7.450003E-02 +1.184011E-03 +3.086286E+00 +1.981929E+00 +1.016918E-01 +2.146313E-03 +4.741530E+00 +4.626144E+00 +4.593581E-02 +4.489950E-04 +1.138847E-01 +2.759746E-03 +2.871515E+00 +1.672027E+00 +8.533571E-02 +1.469319E-03 +4.221958E+00 +3.599412E+00 +3.658547E-02 +2.721530E-04 +9.070318E-02 +1.672788E-03 +2.715053E+00 +1.498990E+00 +8.595576E-02 +1.530999E-03 +4.139891E+00 +3.493193E+00 +3.815373E-02 +3.059779E-04 +9.459122E-02 +1.880692E-03 +2.265823E+00 +1.080584E+00 +6.709890E-02 +9.696052E-04 +3.397886E+00 +2.446009E+00 +2.882467E-02 +1.820798E-04 +7.146249E-02 +1.119153E-03 +2.443336E+00 +1.255468E+00 +8.222450E-02 +1.373726E-03 +3.836077E+00 +3.019588E+00 +3.754438E-02 +2.866712E-04 +9.308051E-02 +1.762024E-03 +2.450506E+00 +1.305510E+00 +7.897423E-02 +1.293368E-03 +3.725732E+00 +2.954033E+00 +3.530661E-02 +2.549027E-04 +8.753260E-02 +1.566759E-03 +2.463151E+00 +1.385265E+00 +8.072623E-02 +1.376083E-03 +3.828575E+00 +3.178469E+00 +3.645911E-02 +2.733269E-04 +9.038989E-02 +1.680003E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -511,86 +511,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.761165E+00 -1.657839E+00 -7.810542E-02 -1.334149E-03 -4.026290E+00 -3.495574E+00 -3.264226E-02 -2.359256E-04 -8.092710E-02 -1.450116E-03 -2.601284E+00 -1.465130E+00 -7.485486E-02 -1.239170E-03 -3.829638E+00 -3.199488E+00 -3.159745E-02 -2.288567E-04 -7.833681E-02 -1.406667E-03 -2.419174E+00 -1.228134E+00 -7.422066E-02 -1.117508E-03 -3.667795E+00 -2.787048E+00 -3.244938E-02 -2.129012E-04 -8.044893E-02 -1.308597E-03 -2.491535E+00 -1.259222E+00 -7.470359E-02 -1.133979E-03 -3.773164E+00 -2.887196E+00 -3.234132E-02 -2.142511E-04 -8.018101E-02 -1.316894E-03 -2.430358E+00 -1.305969E+00 -7.236826E-02 -1.077801E-03 -3.691969E+00 -2.914191E+00 -3.123416E-02 -2.001679E-04 -7.743613E-02 -1.230332E-03 -3.133682E+00 -1.984584E+00 -8.337467E-02 -1.392595E-03 -4.451265E+00 -3.978527E+00 -3.354318E-02 -2.282172E-04 -8.316070E-02 -1.402736E-03 -2.662704E+00 -1.448336E+00 -7.530385E-02 -1.164529E-03 -3.880423E+00 -3.054138E+00 -3.144058E-02 -2.085301E-04 -7.794791E-02 -1.281730E-03 -2.130205E+00 -9.297660E-01 -6.583066E-02 -8.790174E-04 -3.226341E+00 -2.120506E+00 -2.889944E-02 -1.698754E-04 -7.164788E-02 -1.044139E-03 +2.433621E+00 +1.298731E+00 +7.126855E-02 +1.084489E-03 +3.556720E+00 +2.723689E+00 +3.032075E-02 +1.974495E-04 +7.517159E-02 +1.213623E-03 +2.717435E+00 +1.549804E+00 +7.064960E-02 +1.006179E-03 +3.808148E+00 +2.988456E+00 +2.796674E-02 +1.570090E-04 +6.933552E-02 +9.650556E-04 +2.679858E+00 +1.486542E+00 +8.030738E-02 +1.319629E-03 +4.018748E+00 +3.302325E+00 +3.467552E-02 +2.467113E-04 +8.596801E-02 +1.516410E-03 +2.404649E+00 +1.193262E+00 +6.941086E-02 +9.897421E-04 +3.496702E+00 +2.518136E+00 +2.928903E-02 +1.763988E-04 +7.261375E-02 +1.084234E-03 +2.147183E+00 +9.523277E-01 +6.709152E-02 +9.314107E-04 +3.232389E+00 +2.147105E+00 +2.954891E-02 +1.829794E-04 +7.325805E-02 +1.124682E-03 +2.636339E+00 +1.437307E+00 +7.792975E-02 +1.238873E-03 +3.915670E+00 +3.133618E+00 +3.338754E-02 +2.290709E-04 +8.277483E-02 +1.407984E-03 +2.771906E+00 +1.608332E+00 +8.667994E-02 +1.535886E-03 +4.183761E+00 +3.623671E+00 +3.819000E-02 +2.976177E-04 +9.468114E-02 +1.829306E-03 +2.269663E+00 +1.042807E+00 +7.391262E-02 +1.128121E-03 +3.476219E+00 +2.457866E+00 +3.321742E-02 +2.316107E-04 +8.235305E-02 +1.423595E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -681,86 +681,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.599461E+00 -1.458035E+00 -7.809884E-02 -1.305214E-03 -3.901359E+00 -3.258238E+00 -3.379142E-02 -2.432493E-04 -8.377612E-02 -1.495131E-03 -2.748044E+00 -1.595569E+00 -8.641555E-02 -1.591494E-03 -4.133038E+00 -3.610482E+00 -3.815976E-02 -3.125423E-04 -9.460617E-02 -1.921040E-03 -2.827693E+00 -1.681951E+00 -9.214026E-02 -1.815952E-03 -4.343252E+00 -3.980138E+00 -4.139309E-02 -3.763771E-04 -1.026223E-01 -2.313401E-03 -2.848302E+00 -1.709266E+00 -8.618471E-02 -1.555247E-03 -4.293029E+00 -3.893924E+00 -3.742129E-02 -2.947160E-04 -9.277535E-02 -1.811471E-03 -3.107129E+00 -1.985180E+00 -9.390982E-02 -1.778566E-03 -4.705512E+00 -4.502266E+00 -4.077729E-02 -3.361217E-04 -1.010956E-01 -2.065971E-03 -3.176743E+00 -2.076916E+00 -9.535836E-02 -1.843238E-03 -4.782883E+00 -4.655802E+00 -4.125628E-02 -3.453802E-04 -1.022831E-01 -2.122878E-03 -3.158364E+00 -2.020031E+00 -9.047734E-02 -1.652565E-03 -4.641857E+00 -4.326900E+00 -3.809128E-02 -2.983406E-04 -9.443638E-02 -1.833749E-03 -2.750064E+00 -1.515295E+00 -7.550745E-02 -1.157933E-03 -4.022702E+00 -3.249004E+00 -3.106532E-02 -2.005988E-04 -7.701755E-02 -1.232980E-03 +2.364354E+00 +1.160200E+00 +6.641096E-02 +8.852085E-04 +3.392824E+00 +2.329004E+00 +2.757158E-02 +1.580033E-04 +6.835582E-02 +9.711667E-04 +2.524962E+00 +1.341947E+00 +7.077935E-02 +1.024802E-03 +3.665788E+00 +2.770692E+00 +2.942129E-02 +1.757498E-04 +7.294165E-02 +1.080245E-03 +2.400937E+00 +1.186871E+00 +6.912946E-02 +9.715758E-04 +3.511555E+00 +2.501540E+00 +2.915335E-02 +1.747298E-04 +7.227736E-02 +1.073976E-03 +2.114087E+00 +8.990260E-01 +7.043720E-02 +1.000660E-03 +3.333143E+00 +2.235823E+00 +3.206435E-02 +2.087641E-04 +7.949434E-02 +1.283168E-03 +2.130919E+00 +9.474100E-01 +7.007826E-02 +9.974302E-04 +3.308019E+00 +2.242723E+00 +3.167506E-02 +2.044420E-04 +7.852922E-02 +1.256602E-03 +2.315558E+00 +1.159276E+00 +6.983615E-02 +1.023301E-03 +3.433129E+00 +2.511229E+00 +3.018652E-02 +1.890787E-04 +7.483880E-02 +1.162171E-03 +2.729503E+00 +1.530826E+00 +8.250917E-02 +1.373688E-03 +4.093393E+00 +3.424963E+00 +3.577102E-02 +2.569031E-04 +8.868397E-02 +1.579054E-03 +2.303613E+00 +1.077186E+00 +6.880069E-02 +9.471902E-04 +3.366481E+00 +2.277404E+00 +2.953353E-02 +1.751404E-04 +7.321992E-02 +1.076500E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -851,86 +851,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.571690E+00 -1.405222E+00 -6.998368E-02 -1.010782E-03 -3.645470E+00 -2.775516E+00 -2.848750E-02 -1.694087E-04 -7.062657E-02 -1.041270E-03 -3.305320E+00 -2.344814E+00 -9.985328E-02 -2.111130E-03 -4.992554E+00 -5.332303E+00 -4.332586E-02 -3.970490E-04 -1.074140E-01 -2.440461E-03 -2.746305E+00 -1.652960E+00 -9.402849E-02 -1.982149E-03 -4.315533E+00 -4.099428E+00 -4.322622E-02 -4.235148E-04 -1.071670E-01 -2.603132E-03 -2.975372E+00 -1.824531E+00 -9.627993E-02 -1.863677E-03 -4.624041E+00 -4.351874E+00 -4.324350E-02 -3.743520E-04 -1.072099E-01 -2.300953E-03 -2.968077E+00 -1.807761E+00 -9.337061E-02 -1.790184E-03 -4.455263E+00 -4.056544E+00 -4.122926E-02 -3.540009E-04 -1.022161E-01 -2.175865E-03 -3.337278E+00 -2.233238E+00 -1.017763E-01 -2.083989E-03 -4.981699E+00 -4.981154E+00 -4.428100E-02 -3.964318E-04 -1.097820E-01 -2.436667E-03 -3.083638E+00 -1.958766E+00 -9.111141E-02 -1.726305E-03 -4.546082E+00 -4.260420E+00 -3.895365E-02 -3.200133E-04 -9.657439E-02 -1.966961E-03 -2.644733E+00 -1.407404E+00 -8.635553E-02 -1.527997E-03 -4.094426E+00 -3.382785E+00 -3.888081E-02 -3.161126E-04 -9.639381E-02 -1.942985E-03 +2.603449E+00 +1.382310E+00 +7.044110E-02 +1.048854E-03 +3.683939E+00 +2.770822E+00 +2.858207E-02 +1.846566E-04 +7.086104E-02 +1.134991E-03 +2.458052E+00 +1.211571E+00 +7.271340E-02 +1.061659E-03 +3.628656E+00 +2.636961E+00 +3.113049E-02 +1.956492E-04 +7.717911E-02 +1.202557E-03 +2.200265E+00 +9.745445E-01 +7.058665E-02 +1.016272E-03 +3.372745E+00 +2.298238E+00 +3.151097E-02 +2.069373E-04 +7.812240E-02 +1.271940E-03 +2.154223E+00 +9.469347E-01 +6.893166E-02 +1.008453E-03 +3.298063E+00 +2.220933E+00 +3.075978E-02 +2.125545E-04 +7.626005E-02 +1.306466E-03 +2.057172E+00 +8.672290E-01 +6.372568E-02 +8.190278E-04 +3.108552E+00 +1.951197E+00 +2.795535E-02 +1.595890E-04 +6.930727E-02 +9.809131E-04 +2.295642E+00 +1.124640E+00 +6.209725E-02 +8.275041E-04 +3.320499E+00 +2.339919E+00 +2.530846E-02 +1.388437E-04 +6.274506E-02 +8.534023E-04 +2.411516E+00 +1.208668E+00 +7.127755E-02 +1.051206E-03 +3.599456E+00 +2.675604E+00 +3.057946E-02 +1.929511E-04 +7.581300E-02 +1.185973E-03 +2.170599E+00 +9.911770E-01 +6.225331E-02 +8.096954E-04 +3.174400E+00 +2.102442E+00 +2.621745E-02 +1.490631E-04 +6.499865E-02 +9.162161E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -1021,86 +1021,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -3.316522E+00 -2.302876E+00 -8.987919E-02 -1.676655E-03 -4.716109E+00 -4.596488E+00 -3.654407E-02 -2.855910E-04 -9.060053E-02 -1.755384E-03 -3.091331E+00 -2.117646E+00 -8.528196E-02 -1.564336E-03 -4.492326E+00 -4.405850E+00 -3.513814E-02 -2.665883E-04 -8.711494E-02 -1.638584E-03 -2.649730E+00 -1.444519E+00 -8.510948E-02 -1.506094E-03 -4.117595E+00 -3.501517E+00 -3.810553E-02 -3.052819E-04 -9.447172E-02 -1.876414E-03 -2.875773E+00 -1.758531E+00 -9.127582E-02 -1.780419E-03 -4.328266E+00 -3.989840E+00 -4.045386E-02 -3.513542E-04 -1.002937E-01 -2.159598E-03 -3.102792E+00 -1.949906E+00 -9.153879E-02 -1.691646E-03 -4.578530E+00 -4.235906E+00 -3.913672E-02 -3.094156E-04 -9.702826E-02 -1.901822E-03 -3.238743E+00 -2.146355E+00 -8.902551E-02 -1.593536E-03 -4.683916E+00 -4.438028E+00 -3.660342E-02 -2.683183E-04 -9.074766E-02 -1.649218E-03 -3.006635E+00 -1.887385E+00 -8.716712E-02 -1.586834E-03 -4.383965E+00 -4.008888E+00 -3.688262E-02 -2.862513E-04 -9.143987E-02 -1.759443E-03 -2.749904E+00 -1.550561E+00 -9.244273E-02 -1.748082E-03 -4.241273E+00 -3.669144E+00 -4.208622E-02 -3.633241E-04 -1.043407E-01 -2.233170E-03 +2.374348E+00 +1.146696E+00 +6.443426E-02 +8.746572E-04 +3.428224E+00 +2.405912E+00 +2.629872E-02 +1.523681E-04 +6.520012E-02 +9.365302E-04 +2.464893E+00 +1.229701E+00 +7.050588E-02 +1.022217E-03 +3.660115E+00 +2.719644E+00 +2.971255E-02 +1.859463E-04 +7.366374E-02 +1.142918E-03 +2.086598E+00 +8.776819E-01 +6.304625E-02 +8.239181E-04 +3.124875E+00 +1.968039E+00 +2.731211E-02 +1.609716E-04 +6.771252E-02 +9.894116E-04 +2.314873E+00 +1.083111E+00 +5.761556E-02 +6.762901E-04 +3.205662E+00 +2.074555E+00 +2.215216E-02 +1.044406E-04 +5.491992E-02 +6.419437E-04 +2.313273E+00 +1.119912E+00 +6.126598E-02 +7.845228E-04 +3.266340E+00 +2.226545E+00 +2.455456E-02 +1.271226E-04 +6.087598E-02 +7.813588E-04 +2.283282E+00 +1.100474E+00 +6.665553E-02 +9.193170E-04 +3.364565E+00 +2.364829E+00 +2.834631E-02 +1.672037E-04 +7.027653E-02 +1.027717E-03 +2.186072E+00 +9.712029E-01 +7.000241E-02 +1.019717E-03 +3.374290E+00 +2.313941E+00 +3.125803E-02 +2.124583E-04 +7.749531E-02 +1.305874E-03 +2.250156E+00 +1.052173E+00 +6.400255E-02 +8.818690E-04 +3.268015E+00 +2.245559E+00 +2.679559E-02 +1.634377E-04 +6.643197E-02 +1.004569E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -1191,86 +1191,86 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.556952E+00 -1.416458E+00 -8.087467E-02 -1.431667E-03 -3.806748E+00 -3.158182E+00 -3.571995E-02 -2.831832E-04 -8.855737E-02 -1.740585E-03 -2.546384E+00 -1.446967E+00 -8.000281E-02 -1.363415E-03 -3.844739E+00 -3.222133E+00 -3.533522E-02 -2.621956E-04 -8.760354E-02 -1.611584E-03 -2.474418E+00 -1.302438E+00 -7.728372E-02 -1.265047E-03 -3.818517E+00 -3.076074E+00 -3.414744E-02 -2.483215E-04 -8.465877E-02 -1.526307E-03 -2.478272E+00 -1.299434E+00 -7.866430E-02 -1.293627E-03 -3.758685E+00 -2.973817E+00 -3.491240E-02 -2.542033E-04 -8.655529E-02 -1.562460E-03 -2.941937E+00 -1.842278E+00 -9.491310E-02 -1.902183E-03 -4.564730E+00 -4.427111E+00 -4.252722E-02 -3.828960E-04 -1.054340E-01 -2.353469E-03 -2.850564E+00 -1.719152E+00 -8.118192E-02 -1.404765E-03 -4.190858E+00 -3.722715E+00 -3.411695E-02 -2.522778E-04 -8.458318E-02 -1.550625E-03 -2.726097E+00 -1.648250E+00 -7.879124E-02 -1.424955E-03 -3.940978E+00 -3.461222E+00 -3.322044E-02 -2.688615E-04 -8.236053E-02 -1.652557E-03 -2.304822E+00 -1.153775E+00 -7.500010E-02 -1.314077E-03 -3.542935E+00 -2.786864E+00 -3.369367E-02 -2.798215E-04 -8.353377E-02 -1.719922E-03 +2.149693E+00 +9.566763E-01 +6.325375E-02 +8.781459E-04 +3.163781E+00 +2.089993E+00 +2.699802E-02 +1.754950E-04 +6.693384E-02 +1.078679E-03 +2.368262E+00 +1.206255E+00 +6.517665E-02 +9.363557E-04 +3.403305E+00 +2.501730E+00 +2.677251E-02 +1.620249E-04 +6.637476E-02 +9.958857E-04 +2.280249E+00 +1.070455E+00 +7.456831E-02 +1.147671E-03 +3.475673E+00 +2.476848E+00 +3.352210E-02 +2.345869E-04 +8.310842E-02 +1.441888E-03 +2.226928E+00 +1.011037E+00 +7.770236E-02 +1.243393E-03 +3.524398E+00 +2.538453E+00 +3.601410E-02 +2.702997E-04 +8.928663E-02 +1.661396E-03 +2.407411E+00 +1.170093E+00 +7.171133E-02 +1.029235E-03 +3.502484E+00 +2.461827E+00 +3.070568E-02 +1.889715E-04 +7.612591E-02 +1.161513E-03 +2.331359E+00 +1.134448E+00 +6.734317E-02 +9.407363E-04 +3.333577E+00 +2.267345E+00 +2.832929E-02 +1.772720E-04 +7.023433E-02 +1.089602E-03 +2.008838E+00 +8.722179E-01 +5.586199E-02 +7.051708E-04 +2.903675E+00 +1.807255E+00 +2.309145E-02 +1.322637E-04 +5.724863E-02 +8.129585E-04 +2.226061E+00 +1.048269E+00 +6.536859E-02 +9.623671E-04 +3.297988E+00 +2.324450E+00 +2.791406E-02 +1.866630E-04 +6.920490E-02 +1.147324E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -2892,15 +2892,15 @@ tally 1: 0.000000E+00 0.000000E+00 tally 2: -4.283244E+01 -3.698890E+02 -4.285612E+01 -3.702981E+02 -6.926001E+00 -9.669342E+00 -6.926497E+00 -9.670727E+00 -1.223563E+02 -3.025594E+03 -1.223563E+02 -3.025594E+03 +4.076711E+01 +3.341374E+02 +4.077838E+01 +3.343221E+02 +6.274781E+00 +7.937573E+00 +6.275007E+00 +7.938146E+00 +1.122968E+02 +2.557771E+03 +1.122968E+02 +2.557771E+03