diff --git a/CMakeLists.txt b/CMakeLists.txt index a7e8da55a7..100e148688 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -378,6 +378,7 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/tally.F90 src/tallies/tally_filter.F90 src/tallies/tally_filter_header.F90 + src/tallies/tally_filter_energy.F90 src/tallies/tally_filter_mesh.F90 src/tallies/tally_header.F90 src/tallies/tally_initialize.F90 diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 55862cd522..2a9f1da306 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -254,7 +254,6 @@ contains use tally_header, only: TallyObject use tally_filter_header use tally_filter - use tally_filter_mesh, only: MeshFilter use tally_initialize, only: add_tallies use xml_interface diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6923b39349..911d2fe554 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -30,7 +30,6 @@ module input_xml use tally_header, only: TallyObject use tally_filter_header, only: TallyFilterContainer use tally_filter - use tally_filter_mesh, only: MeshFilter use tally_initialize, only: add_tallies use xml_interface diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 404a1a4877..2472143651 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -15,7 +15,6 @@ module tally use particle_header, only: LocalCoord, Particle use string, only: to_str use tally_filter - use tally_filter_mesh, only: MeshFilter implicit none diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index a718202b78..4b1e5acb94 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -11,7 +11,11 @@ module tally_filter use particle_header, only: Particle use string, only: to_str use tally_filter_header, only: TallyFilter, TallyFilterContainer, & - TallyFilterMatch + TallyFilterMatch + + ! Inherit other filters + use tally_filter_energy + use tally_filter_mesh implicit none @@ -106,38 +110,6 @@ module tally_filter procedure :: initialize => initialize_surface end type SurfaceFilter -!=============================================================================== -! ENERGYFILTER bins the incident neutron energy. -!=============================================================================== - type, extends(TallyFilter) :: EnergyFilter - real(8), allocatable :: bins(:) - - ! True if transport group number can be used directly to get bin number - logical :: matches_transport_groups = .false. - - contains - procedure :: get_all_bins => get_all_bins_energy - procedure :: to_statepoint => to_statepoint_energy - procedure :: text_label => text_label_energy - end type EnergyFilter - -!=============================================================================== -! ENERGYOUTFILTER bins the outgoing neutron energy. Only scattering events use -! the get_all_bins functionality. Nu-fission tallies manually iterate over the -! filter bins. -!=============================================================================== - type, extends(TallyFilter) :: EnergyoutFilter - real(8), allocatable :: bins(:) - - ! True if transport group number can be used directly to get bin number - logical :: matches_transport_groups = .false. - - contains - procedure :: get_all_bins => get_all_bins_energyout - procedure :: to_statepoint => to_statepoint_energyout - procedure :: text_label => text_label_energyout - end type EnergyoutFilter - !=============================================================================== ! DELAYEDGROUPFILTER bins outgoing fission neutrons in their delayed groups. ! The get_all_bins functionality is not actually used. The bins are manually @@ -653,118 +625,6 @@ contains label = "Surface " // to_str(surfaces(this % surfaces(bin)) % obj % id) end function text_label_surface -!=============================================================================== -! EnergyFilter methods -!=============================================================================== - subroutine get_all_bins_energy(this, p, estimator, match) - class(EnergyFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: E - - n = this % n_bins - - if (p % g /= NONE .and. this % matches_transport_groups) then - if (estimator == ESTIMATOR_TRACKLENGTH) then - call match % bins % push_back(num_energy_groups - p % g + 1) - call match % weights % push_back(ONE) - else - call match % bins % push_back(num_energy_groups - p % last_g + 1) - call match % weights % push_back(ONE) - end if - - else - ! Pre-collision energy of particle - E = p % last_E - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, E) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - end if - end subroutine get_all_bins_energy - - subroutine to_statepoint_energy(this, filter_group) - class(EnergyFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "energy") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_energy - - function text_label_energy(this, bin) result(label) - class(EnergyFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Incoming Energy [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_energy - -!=============================================================================== -! EnergyoutFilter methods -!=============================================================================== - subroutine get_all_bins_energyout(this, p, estimator, match) - class(EnergyoutFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - - n = this % n_bins - - if (p % g /= NONE .and. this % matches_transport_groups) then - ! Tallies are ordered in increasing groups, group indices - ! however are the opposite, so switch - call match % bins % push_back(num_energy_groups - p % g + 1) - call match % weights % push_back(ONE) - - else - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, p % E) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - end if - end subroutine get_all_bins_energyout - - subroutine to_statepoint_energyout(this, filter_group) - class(EnergyoutFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "energyout") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_energyout - - function text_label_energyout(this, bin) result(label) - class(EnergyoutFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Outgoing Energy [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_energyout - !=============================================================================== ! DelayedGroupFilter methods !=============================================================================== diff --git a/src/tallies/tally_filter_energy.F90 b/src/tallies/tally_filter_energy.F90 new file mode 100644 index 0000000000..29a203d17c --- /dev/null +++ b/src/tallies/tally_filter_energy.F90 @@ -0,0 +1,160 @@ +module tally_filter_energy + + use hdf5, only: HID_T + + use algorithm, only: binary_search + use constants + use hdf5_interface + use mgxs_header, only: num_energy_groups + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header, only: TallyFilter, TallyFilterMatch + + implicit none + private + +!=============================================================================== +! ENERGYFILTER bins the incident neutron energy. +!=============================================================================== + + type, public, extends(TallyFilter) :: EnergyFilter + real(8), allocatable :: bins(:) + + ! True if transport group number can be used directly to get bin number + logical :: matches_transport_groups = .false. + contains + procedure :: get_all_bins => get_all_bins_energy + procedure :: to_statepoint => to_statepoint_energy + procedure :: text_label => text_label_energy + end type EnergyFilter + +!=============================================================================== +! ENERGYOUTFILTER bins the outgoing neutron energy. Only scattering events use +! the get_all_bins functionality. Nu-fission tallies manually iterate over the +! filter bins. +!=============================================================================== + + type, public, extends(EnergyFilter) :: EnergyoutFilter + contains + procedure :: get_all_bins => get_all_bins_energyout + procedure :: to_statepoint => to_statepoint_energyout + procedure :: text_label => text_label_energyout + end type EnergyoutFilter + +contains + +!=============================================================================== +! EnergyFilter methods +!=============================================================================== + + subroutine get_all_bins_energy(this, p, estimator, match) + class(EnergyFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: n + integer :: bin + real(8) :: E + + n = this % n_bins + + if (p % g /= NONE .and. this % matches_transport_groups) then + if (estimator == ESTIMATOR_TRACKLENGTH) then + call match % bins % push_back(num_energy_groups - p % g + 1) + call match % weights % push_back(ONE) + else + call match % bins % push_back(num_energy_groups - p % last_g + 1) + call match % weights % push_back(ONE) + end if + + else + ! Pre-collision energy of particle + E = p % last_E + + ! Search to find incoming energy bin. + bin = binary_search(this % bins, n + 1, E) + if (bin /= NO_BIN_FOUND) then + call match % bins % push_back(bin) + call match % weights % push_back(ONE) + end if + end if + end subroutine get_all_bins_energy + + subroutine to_statepoint_energy(this, filter_group) + class(EnergyFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "energy") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % bins) + end subroutine to_statepoint_energy + + function text_label_energy(this, bin) result(label) + class(EnergyFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + real(8) :: E0, E1 + + E0 = this % bins(bin) + E1 = this % bins(bin + 1) + label = "Incoming Energy [" // trim(to_str(E0)) // ", " & + // trim(to_str(E1)) // ")" + end function text_label_energy + +!=============================================================================== +! EnergyoutFilter methods +!=============================================================================== + + subroutine get_all_bins_energyout(this, p, estimator, match) + class(EnergyoutFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: n + integer :: bin + + n = this % n_bins + + if (p % g /= NONE .and. this % matches_transport_groups) then + ! Tallies are ordered in increasing groups, group indices + ! however are the opposite, so switch + call match % bins % push_back(num_energy_groups - p % g + 1) + call match % weights % push_back(ONE) + + else + + ! Search to find incoming energy bin. + bin = binary_search(this % bins, n + 1, p % E) + if (bin /= NO_BIN_FOUND) then + call match % bins % push_back(bin) + call match % weights % push_back(ONE) + end if + end if + end subroutine get_all_bins_energyout + + subroutine to_statepoint_energyout(this, filter_group) + class(EnergyoutFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "energyout") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % bins) + end subroutine to_statepoint_energyout + + function text_label_energyout(this, bin) result(label) + class(EnergyoutFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + real(8) :: E0, E1 + + E0 = this % bins(bin) + E1 = this % bins(bin + 1) + label = "Outgoing Energy [" // trim(to_str(E0)) // ", " & + // trim(to_str(E1)) // ")" + end function text_label_energyout + +end module tally_filter_energy diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 5c8c3ba095..f1c8615f57 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -11,9 +11,7 @@ module tally_filter_mesh use tally_filter_header, only: TallyFilter, TallyFilterMatch implicit none - private - public :: MeshFilter !=============================================================================== ! MESHFILTER indexes the location of particle events to a regular mesh. For @@ -21,7 +19,7 @@ module tally_filter_mesh ! will correspond to the fraction of the track length that lies in that bin. !=============================================================================== - type, extends(TallyFilter) :: MeshFilter + type, public, extends(TallyFilter) :: MeshFilter type(RegularMesh), pointer :: mesh => null() contains procedure :: get_all_bins => get_all_bins_mesh