crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Fixed merge conflicts with develop

Browse source
This commit is contained in:
Will Boyd 2016-05-12 21:27:09 -04:00
commit 223e22d9ae
27 changed files with 3503 additions and 2998 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/ace.F90
View file

@ -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

8
src/cross_section.F90
View file

@ -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

6
src/energy_grid.F90
View file

@ -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

8
src/global.F90
View file

@ -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

57
src/input_xml.F90
View file

@ -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

860
src/macroxs_header.F90
View file

@ -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

27
src/math.F90
View file

@ -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

58
src/mgxs_data.F90
View file

@ -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

1929
src/mgxs_header.F90 Normal file
View file

File diff suppressed because it is too large Load diff

898
src/nuclide_header.F90
View file

@ -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

19
src/output.F90
View file

@ -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

2
src/particle_header.F90
View file

@ -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

24
src/physics.F90
View file

@ -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

15
src/physics_mg.F90
View file

@ -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

697
src/scattdata_header.F90
View file

@ -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

6
src/state_point.F90
View file

@ -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.

6
src/summary.F90
View file

@ -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)

534
src/tally.F90
View file

@ -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

2
src/tracking.F90
View file

@ -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

2
tests/test_mg_basic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.045320E+00 5.851680E-02
1.033731E+00 4.974463E-02

2
tests/test_mg_max_order/inputs_true.dat
View file

@ -1 +1 @@
7f7465abaf559b3ef56cb6b0f28050c12f392f55db33dc5d2cefc14b92beb2c9068834c05273e51323d3516643e8a385e4c177a7a471678c961808d19055a30f
60a35864ad71646309d7f1687ba0826d4d53a5b2e8babf73614362645205484bad3c0e7bf605ec0b11cadf58474b2e3d0a97bf2d9297f9118682c37ff0269afd

2
tests/test_mg_max_order/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.083030E+00 1.855038E-02
1.055274E+00 1.715904E-02

3
tests/test_mg_max_order/test_mg_max_order.py
View file

@ -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)

2
tests/test_mg_nuclide/inputs_true.dat
View file

@ -1 +1 @@
825dee3ca35d48788f1a4d5364789bbd83b36e33af9a990da758dd73c3bfcbee14bce2a41e6c80e0147f45575e59078653c8dfa8590cd361c09f19c26dc8c88e
0efba3dd7882fdd38756d0a8f01ff00d7a1abdaab6430b3f090f3339e552448453bbb733852b6bd6ff09608d923c282f168320f942fc2eb3a45610873c588734

2
tests/test_mg_nuclide/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.380785E-01 5.556526E-03
1.033731E+00 4.974463E-02

8
tests/test_mg_nuclide/test_mg_nuclide.py
View file

@ -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.

1306
tests/test_mg_tallies/results_true.dat
View file

File diff suppressed because it is too large Load diff