diff --git a/CMakeLists.txt b/CMakeLists.txt index 86bc20f59..168663110 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -119,8 +119,8 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU) list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2) if(debug) list(REMOVE_ITEM f90flags -O2) - list(APPEND f90flags -g -Wall -pedantic -fbounds-check - -ffpe-trap=invalid,overflow,underflow) + list(APPEND f90flags -g -Wall -Wno-unused-dummy-argument -pedantic + -fbounds-check -ffpe-trap=invalid,overflow,underflow) list(APPEND ldflags -g) endif() if(profile) @@ -379,10 +379,20 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/tally.F90 src/tallies/tally_filter.F90 src/tallies/tally_filter_header.F90 + src/tallies/tally_filter_azimuthal.F90 + src/tallies/tally_filter_cell.F90 + src/tallies/tally_filter_cellborn.F90 + src/tallies/tally_filter_cellfrom.F90 + src/tallies/tally_filter_delayedgroup.F90 src/tallies/tally_filter_distribcell.F90 src/tallies/tally_filter_energy.F90 + src/tallies/tally_filter_energyfunc.F90 src/tallies/tally_filter_material.F90 src/tallies/tally_filter_mesh.F90 + src/tallies/tally_filter_mu.F90 + src/tallies/tally_filter_polar.F90 + src/tallies/tally_filter_surface.F90 + src/tallies/tally_filter_universe.F90 src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 diff --git a/src/settings.F90 b/src/settings.F90 new file mode 100644 index 000000000..d1d41a75c --- /dev/null +++ b/src/settings.F90 @@ -0,0 +1,141 @@ +module settings + + use constants + use set_header, only: SetInt + use source_header + + implicit none + + ! ============================================================================ + ! ENERGY TREATMENT RELATED VARIABLES + logical :: run_CE = .true. ! Run in CE mode? + + ! ============================================================================ + ! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES + + ! Unreoslved resonance probablity tables + logical :: urr_ptables_on = .true. + + ! Default temperature and method for choosing temperatures + integer :: temperature_method = TEMPERATURE_NEAREST + logical :: temperature_multipole = .false. + real(8) :: temperature_tolerance = 10.0_8 + real(8) :: temperature_default = 293.6_8 + real(8) :: temperature_range(2) = [ZERO, ZERO] + + integer :: n_log_bins ! number of bins for logarithmic grid + + ! ============================================================================ + ! MULTI-GROUP CROSS SECTION RELATED VARIABLES + + ! Maximum Data Order + integer :: max_order + + ! Whether or not to convert Legendres to tabulars + logical :: legendre_to_tabular = .true. + + ! Number of points to use in the Legendre to tabular conversion + integer :: legendre_to_tabular_points = 33 + + ! Assume all tallies are spatially distinct + logical :: assume_separate = .false. + + ! Use confidence intervals for results instead of standard deviations + logical :: confidence_intervals = .false. + + ! ============================================================================ + ! SIMULATION VARIABLES + + integer(8) :: n_particles = 0 ! # of particles per generation + integer :: n_batches ! # of batches + integer :: n_inactive ! # of inactive batches + integer :: n_active ! # of active batches + integer :: gen_per_batch = 1 ! # of generations per batch + + integer :: n_max_batches ! max # of batches + integer :: n_batch_interval = 1 ! batch interval for triggers + logical :: pred_batches = .false. ! predict batches for triggers + logical :: trigger_on = .false. ! flag for turning triggers on/off + + logical :: entropy_on = .false. + integer :: index_entropy_mesh = -1 + + logical :: ufs = .false. + integer :: index_ufs_mesh = -1 + + ! Write source at end of simulation + logical :: source_separate = .false. + logical :: source_write = .true. + logical :: source_latest = .false. + + ! Variance reduction settins + logical :: survival_biasing = .false. + real(8) :: weight_cutoff = 0.25_8 + real(8) :: energy_cutoff = ZERO + real(8) :: weight_survive = ONE + + ! Mode to run in (fixed source, eigenvalue, plotting, etc) + integer :: run_mode = NONE + + ! Restart run + logical :: restart_run = .false. + + ! The verbosity controls how much information will be printed to the screen + ! and in logs + integer :: verbosity = 7 + + logical :: check_overlaps = .false. + + ! Trace for single particle + integer :: trace_batch + integer :: trace_gen + integer(8) :: trace_particle + + ! Particle tracks + logical :: write_all_tracks = .false. + integer, allocatable :: track_identifiers(:,:) + + ! Particle restart run + logical :: particle_restart_run = .false. + + ! Write out initial source + logical :: write_initial_source = .false. + + ! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE + logical :: create_fission_neutrons = .true. + + ! Information about state points to be written + integer :: n_state_points = 0 + type(SetInt) :: statepoint_batch + + ! Information about source points to be written + integer :: n_source_points = 0 + type(SetInt) :: sourcepoint_batch + + character(MAX_FILE_LEN) :: path_input ! Path to input file + character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml + character(MAX_FILE_LEN) :: path_multipole ! Path to wmp library + character(MAX_FILE_LEN) :: path_source = '' ! Path to binary source + character(MAX_FILE_LEN) :: path_state_point ! Path to binary state point + character(MAX_FILE_LEN) :: path_source_point ! Path to binary source point + character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart + character(MAX_FILE_LEN) :: path_output = '' ! Path to output directory + + ! Various output options + logical :: output_summary = .true. + logical :: output_tallies = .true. + + ! Resonance scattering settings + logical :: res_scat_on = .false. ! is resonance scattering treated? + integer :: res_scat_method = RES_SCAT_ARES ! resonance scattering method + real(8) :: res_scat_energy_min = 0.01_8 + real(8) :: res_scat_energy_max = 1000.0_8 + character(10), allocatable :: res_scat_nuclides(:) + + ! Is CMFD active + logical :: cmfd_run = .false. + + ! No reduction at end of batch + logical :: reduce_tallies = .true. + +end module settings diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index dc3aa36b1..3b96cade3 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -1,704 +1,33 @@ module tally_filter + use, intrinsic :: ISO_C_BINDING + use hdf5, only: HID_T - use algorithm, only: binary_search - use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION, ERROR_REAL - use dict_header, only: DictIntInt use error - use geometry_header - use hdf5_interface - use particle_header, only: Particle - use surface_header - use string, only: to_str, to_f_string + use string, only: to_f_string use tally_filter_header ! Inherit other filters + use tally_filter_azimuthal + use tally_filter_cell + use tally_filter_cellborn + use tally_filter_cellfrom + use tally_filter_delayedgroup use tally_filter_distribcell use tally_filter_energy + use tally_filter_energyfunc use tally_filter_material use tally_filter_mesh + use tally_filter_mu + use tally_filter_polar + use tally_filter_surface + use tally_filter_universe implicit none -!=============================================================================== -! UNIVERSEFILTER specifies which geometric universes tally events reside in. -!=============================================================================== - type, extends(TallyFilter) :: UniverseFilter - integer, allocatable :: universes(:) - type(DictIntInt) :: map - contains - procedure :: get_all_bins => get_all_bins_universe - procedure :: to_statepoint => to_statepoint_universe - procedure :: text_label => text_label_universe - procedure :: initialize => initialize_universe - end type UniverseFilter - -!=============================================================================== -! CELLFILTER specifies which geometric cells tally events reside in. -!=============================================================================== - type, extends(TallyFilter) :: CellFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map - contains - procedure :: get_all_bins => get_all_bins_cell - procedure :: to_statepoint => to_statepoint_cell - procedure :: text_label => text_label_cell - procedure :: initialize => initialize_cell - end type CellFilter - -!=============================================================================== -! CELLFROMFILTER specifies which geometric cells particles exit when crossing a -! surface. -!=============================================================================== - type, extends(CellFilter) :: CellFromFilter - contains - procedure :: get_all_bins => get_all_bins_cell_from - procedure :: to_statepoint => to_statepoint_cell_from - procedure :: text_label => text_label_cell_from - end type CellFromFilter - -!=============================================================================== -! CELLBORNFILTER specifies which cell the particle was born in. -!=============================================================================== - type, extends(TallyFilter) :: CellbornFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map - contains - procedure :: get_all_bins => get_all_bins_cellborn - procedure :: to_statepoint => to_statepoint_cellborn - procedure :: text_label => text_label_cellborn - procedure :: initialize => initialize_cellborn - end type CellbornFilter - -!=============================================================================== -! SURFACEFILTER specifies which surface particles are crossing -!=============================================================================== - type, extends(TallyFilter) :: SurfaceFilter - integer, allocatable :: surfaces(:) - - ! True if this filter is used for surface currents - logical :: current = .false. - contains - procedure :: get_all_bins => get_all_bins_surface - procedure :: to_statepoint => to_statepoint_surface - procedure :: text_label => text_label_surface - procedure :: initialize => initialize_surface - end type SurfaceFilter - -!=============================================================================== -! DELAYEDGROUPFILTER bins outgoing fission neutrons in their delayed groups. -! The get_all_bins functionality is not actually used. The bins are manually -! iterated over in the scoring subroutines. -!=============================================================================== - type, extends(TallyFilter) :: DelayedGroupFilter - integer, allocatable :: groups(:) - contains - procedure :: get_all_bins => get_all_bins_dg - procedure :: to_statepoint => to_statepoint_dg - procedure :: text_label => text_label_dg - end type DelayedGroupFilter - -!=============================================================================== -! MUFILTER bins the incoming-outgoing direction cosine. This is only used for -! scatter reactions. -!=============================================================================== - type, extends(TallyFilter) :: MuFilter - real(8), allocatable :: bins(:) - contains - procedure :: get_all_bins => get_all_bins_mu - procedure :: to_statepoint => to_statepoint_mu - procedure :: text_label => text_label_mu - end type MuFilter - -!=============================================================================== -! POLARFILTER bins the incident neutron polar angle (relative to the global -! z-axis). -!=============================================================================== - type, extends(TallyFilter) :: PolarFilter - real(8), allocatable :: bins(:) - contains - procedure :: get_all_bins => get_all_bins_polar - procedure :: to_statepoint => to_statepoint_polar - procedure :: text_label => text_label_polar - end type PolarFilter - -!=============================================================================== -! AZIMUTHALFILTER bins the incident neutron azimuthal angle (relative to the -! global xy-plane). -!=============================================================================== - type, extends(TallyFilter) :: AzimuthalFilter - real(8), allocatable :: bins(:) - contains - procedure :: get_all_bins => get_all_bins_azimuthal - procedure :: to_statepoint => to_statepoint_azimuthal - procedure :: text_label => text_label_azimuthal - end type AzimuthalFilter - -!=============================================================================== -! EnergyFunctionFilter multiplies tally scores by an arbitrary function of -! incident energy described by a piecewise linear-linear interpolation. -!=============================================================================== - type, extends(TallyFilter) :: EnergyFunctionFilter - real(8), allocatable :: energy(:) - real(8), allocatable :: y(:) - - contains - procedure :: get_all_bins => get_all_bins_energyfunction - procedure :: to_statepoint => to_statepoint_energyfunction - procedure :: text_label => text_label_energyfunction - end type EnergyFunctionFilter - contains -!=============================================================================== -! METHODS: for a description of these methods, see their counterparts bound to -! the abstract TallyFilter class. -!=============================================================================== - -!=============================================================================== -! UniverseFilter methods -!=============================================================================== - subroutine get_all_bins_universe(this, p, estimator, match) - class(UniverseFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - - ! Iterate over coordinate levels to see which universes match - do i = 1, p % n_coord - if (this % map % has_key(p % coord(i) % universe)) then - call match % bins % push_back(this % map % get_key(p % coord(i) & - % universe)) - call match % weights % push_back(ONE) - end if - end do - - end subroutine get_all_bins_universe - - subroutine to_statepoint_universe(this, filter_group) - class(UniverseFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: universe_ids(:) - - call write_dataset(filter_group, "type", "universe") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(universe_ids(size(this % universes))) - do i = 1, size(this % universes) - universe_ids(i) = universes(this % universes(i)) % id - end do - call write_dataset(filter_group, "bins", universe_ids) - end subroutine to_statepoint_universe - - subroutine initialize_universe(this) - class(UniverseFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % universes(i) - if (universe_dict % has_key(id)) then - this % universes(i) = universe_dict % get_key(id) - else - call fatal_error("Could not find universe " // trim(to_str(id)) & - &// " specified on a tally filter.") - end if - end do - - ! Generate mapping from universe indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % universes(i), i) - end do - end subroutine initialize_universe - - function text_label_universe(this, bin) result(label) - class(UniverseFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Universe " // to_str(universes(this % universes(bin)) % id) - end function text_label_universe - -!=============================================================================== -! CellFilter methods -!=============================================================================== - subroutine get_all_bins_cell(this, p, estimator, match) - class(CellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - - ! Iterate over coordinate levels to see with cells match - do i = 1, p % n_coord - if (this % map % has_key(p % coord(i) % cell)) then - call match % bins % push_back(this % map % get_key(p % coord(i) % cell)) - call match % weights % push_back(ONE) - end if - end do - - end subroutine get_all_bins_cell - - subroutine to_statepoint_cell(this, filter_group) - class(CellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cell") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell - - subroutine initialize_cell(this) - class(CellFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) - if (cell_dict % has_key(id)) then - this % cells(i) = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % cells(i), i) - end do - end subroutine initialize_cell - - function text_label_cell(this, bin) result(label) - class(CellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell " // to_str(cells(this % cells(bin)) % id) - end function text_label_cell - -!=============================================================================== -! CellFromFilter methods -!=============================================================================== - subroutine get_all_bins_cell_from(this, p, estimator, match) - class(CellFromFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - - ! Starting one coordinate level deeper, find the next bin. - do i = 1, p % last_n_coord - if (this % map % has_key(p % last_cell(i))) then - call match % bins % push_back(this % map % get_key(p % last_cell(i))) - call match % weights % push_back(ONE) - exit - end if - end do - - end subroutine get_all_bins_cell_from - - subroutine to_statepoint_cell_from(this, filter_group) - class(CellFromFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cellfrom") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell_from - - function text_label_cell_from(this, bin) result(label) - class(CellFromFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell from " // to_str(cells(this % cells(bin)) % id) - end function text_label_cell_from - -!=============================================================================== -! CellbornFilter methods -!=============================================================================== - subroutine get_all_bins_cellborn(this, p, estimator, match) - class(CellbornFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - if (this % map % has_key(p % cell_born)) then - call match % bins % push_back(this % map % get_key(p % cell_born)) - call match % weights % push_back(ONE) - end if - - end subroutine get_all_bins_cellborn - - subroutine to_statepoint_cellborn(this, filter_group) - class(CellbornFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cellborn") - call write_dataset(filter_group, "n_bins", this % n_bins) - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cellborn - - subroutine initialize_cellborn(this) - class(CellbornFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) - if (cell_dict % has_key(id)) then - this % cells(i) = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % cells(i), i) - end do - end subroutine initialize_cellborn - - function text_label_cellborn(this, bin) result(label) - class(CellbornFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Birth Cell " // to_str(cells(this % cells(bin)) % id) - end function text_label_cellborn - -!=============================================================================== -! SurfaceFilter methods -!=============================================================================== - subroutine get_all_bins_surface(this, p, estimator, match) - class(SurfaceFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - - do i = 1, this % n_bins - if (abs(p % surface) == this % surfaces(i)) then - call match % bins % push_back(i) - if (p % surface < 0) then - call match % weights % push_back(-ONE) - else - call match % weights % push_back(ONE) - end if - exit - end if - end do - - end subroutine get_all_bins_surface - - subroutine to_statepoint_surface(this, filter_group) - class(SurfaceFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "surface") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % surfaces) - end subroutine to_statepoint_surface - - subroutine initialize_surface(this) - class(SurfaceFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % surfaces(i) - if (surface_dict % has_key(id)) then - this % surfaces(i) = surface_dict % get_key(id) - else - call fatal_error("Could not find surface " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - end subroutine initialize_surface - - function text_label_surface(this, bin) result(label) - class(SurfaceFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Surface " // to_str(surfaces(this % surfaces(bin)) % obj % id) - end function text_label_surface - -!=============================================================================== -! DelayedGroupFilter methods -!=============================================================================== - subroutine get_all_bins_dg(this, p, estimator, match) - class(DelayedGroupFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - call match % bins % push_back(1) - call match % weights % push_back(ONE) - end subroutine get_all_bins_dg - - subroutine to_statepoint_dg(this, filter_group) - class(DelayedGroupFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "delayedgroup") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % groups) - end subroutine to_statepoint_dg - - function text_label_dg(this, bin) result(label) - class(DelayedGroupFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Delayed Group " // to_str(this % groups(bin)) - end function text_label_dg - -!=============================================================================== -! MuFilter methods -!=============================================================================== - subroutine get_all_bins_mu(this, p, estimator, match) - class(MuFilter), 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 - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, p % mu) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - end subroutine get_all_bins_mu - - subroutine to_statepoint_mu(this, filter_group) - class(MuFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "mu") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_mu - - function text_label_mu(this, bin) result(label) - class(MuFilter), 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 = "Change-in-Angle [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_mu - -!=============================================================================== -! PolarFilter methods -!=============================================================================== - subroutine get_all_bins_polar(this, p, estimator, match) - class(PolarFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: theta - - n = this % n_bins - - ! Make sure the correct direction vector is used. - if (estimator == ESTIMATOR_TRACKLENGTH) then - theta = acos(p % coord(1) % uvw(3)) - else - theta = acos(p % last_uvw(3)) - end if - - ! Search to find polar angle bin. - bin = binary_search(this % bins, n + 1, theta) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - end subroutine get_all_bins_polar - - subroutine to_statepoint_polar(this, filter_group) - class(PolarFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "polar") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_polar - - function text_label_polar(this, bin) result(label) - class(PolarFilter), 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 = "Polar Angle [" // trim(to_str(E0)) // ", " // trim(to_str(E1)) & - // ")" - end function text_label_polar - -!=============================================================================== -! AzimuthalFilter methods -!=============================================================================== - subroutine get_all_bins_azimuthal(this, p, estimator, match) - class(AzimuthalFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: phi - - n = this % n_bins - - ! Make sure the correct direction vector is used. - if (estimator == ESTIMATOR_TRACKLENGTH) then - phi = atan2(p % coord(1) % uvw(2), p % coord(1) % uvw(1)) - else - phi = atan2(p % last_uvw(2), p % last_uvw(1)) - end if - - ! Search to find azimuthal angle bin. - bin = binary_search(this % bins, n + 1, phi) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - - end subroutine get_all_bins_azimuthal - - subroutine to_statepoint_azimuthal(this, filter_group) - class(AzimuthalFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "azimuthal") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_azimuthal - - function text_label_azimuthal(this, bin) result(label) - class(AzimuthalFilter), 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 = "Azimuthal Angle [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_azimuthal - -!=============================================================================== -! EnergyFunctionFilter methods -!=============================================================================== - subroutine get_all_bins_energyfunction(this, p, estimator, match) - class(EnergyFunctionFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n, indx - real(8) :: E, f, weight - - select type(this) - type is (EnergyFunctionFilter) - n = size(this % energy) - - ! Get pre-collision energy of particle - E = p % last_E - - ! Search to find incoming energy bin. - indx = binary_search(this % energy, n, E) - - ! Compute an interpolation factor between nearest bins. - f = (E - this % energy(indx)) & - / (this % energy(indx+1) - this % energy(indx)) - - ! Interpolate on the lin-lin grid. - call match % bins % push_back(1) - weight = (ONE - f) * this % y(indx) + f * this % y(indx+1) - call match % weights % push_back(weight) - end select - end subroutine get_all_bins_energyfunction - - subroutine to_statepoint_energyfunction(this, filter_group) - class(EnergyFunctionFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - select type(this) - type is (EnergyFunctionFilter) - call write_dataset(filter_group, "type", "energyfunction") - call write_dataset(filter_group, "energy", this % energy) - call write_dataset(filter_group, "y", this % y) - end select - end subroutine to_statepoint_energyfunction - - function text_label_energyfunction(this, bin) result(label) - class(EnergyFunctionFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - select type(this) - type is (EnergyFunctionFilter) - write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") & - "Energy Function f([", this % energy(1), ", ..., ", & - this % energy(size(this % energy)), "]) = [", this % y(1), & - ", ..., ", this % y(size(this % y)), "]" - end select - end function text_label_energyfunction - !=============================================================================== ! C API FUNCTIONS !=============================================================================== diff --git a/src/tallies/tally_filter_azimuthal.F90 b/src/tallies/tally_filter_azimuthal.F90 new file mode 100644 index 000000000..a6c92fb8e --- /dev/null +++ b/src/tallies/tally_filter_azimuthal.F90 @@ -0,0 +1,83 @@ +module tally_filter_azimuthal + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use algorithm, only: binary_search + use constants + use error, only: fatal_error + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! AZIMUTHALFILTER bins the incident neutron azimuthal angle (relative to the +! global xy-plane). +!=============================================================================== + + type, public, extends(TallyFilter) :: AzimuthalFilter + real(8), allocatable :: bins(:) + contains + procedure :: get_all_bins => get_all_bins_azimuthal + procedure :: to_statepoint => to_statepoint_azimuthal + procedure :: text_label => text_label_azimuthal + end type AzimuthalFilter + +contains + + subroutine get_all_bins_azimuthal(this, p, estimator, match) + class(AzimuthalFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: n + integer :: bin + real(8) :: phi + + n = this % n_bins + + ! Make sure the correct direction vector is used. + if (estimator == ESTIMATOR_TRACKLENGTH) then + phi = atan2(p % coord(1) % uvw(2), p % coord(1) % uvw(1)) + else + phi = atan2(p % last_uvw(2), p % last_uvw(1)) + end if + + ! Search to find azimuthal angle bin. + bin = binary_search(this % bins, n + 1, phi) + if (bin /= NO_BIN_FOUND) then + call match % bins % push_back(bin) + call match % weights % push_back(ONE) + end if + + end subroutine get_all_bins_azimuthal + + subroutine to_statepoint_azimuthal(this, filter_group) + class(AzimuthalFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "azimuthal") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % bins) + end subroutine to_statepoint_azimuthal + + function text_label_azimuthal(this, bin) result(label) + class(AzimuthalFilter), 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 = "Azimuthal Angle [" // trim(to_str(E0)) // ", " & + // trim(to_str(E1)) // ")" + end function text_label_azimuthal + +end module tally_filter_azimuthal diff --git a/src/tallies/tally_filter_cell.F90 b/src/tallies/tally_filter_cell.F90 new file mode 100644 index 000000000..2171e05e9 --- /dev/null +++ b/src/tallies/tally_filter_cell.F90 @@ -0,0 +1,99 @@ +module tally_filter_cell + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error + use hdf5_interface + use geometry_header + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! CELLFILTER specifies which geometric cells tally events reside in. +!=============================================================================== + + type, public, extends(TallyFilter) :: CellFilter + integer, allocatable :: cells(:) + type(DictIntInt) :: map + contains + procedure :: get_all_bins => get_all_bins_cell + procedure :: to_statepoint => to_statepoint_cell + procedure :: text_label => text_label_cell + procedure :: initialize => initialize_cell + end type CellFilter + +contains + + subroutine get_all_bins_cell(this, p, estimator, match) + class(CellFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + + ! Iterate over coordinate levels to see with cells match + do i = 1, p % n_coord + if (this % map % has_key(p % coord(i) % cell)) then + call match % bins % push_back(this % map % get_key(p % coord(i) % cell)) + call match % weights % push_back(ONE) + end if + end do + + end subroutine get_all_bins_cell + + subroutine to_statepoint_cell(this, filter_group) + class(CellFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: cell_ids(:) + + call write_dataset(filter_group, "type", "cell") + call write_dataset(filter_group, "n_bins", this % n_bins) + + allocate(cell_ids(size(this % cells))) + do i = 1, size(this % cells) + cell_ids(i) = cells(this % cells(i)) % id + end do + call write_dataset(filter_group, "bins", cell_ids) + end subroutine to_statepoint_cell + + subroutine initialize_cell(this) + class(CellFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % cells(i) + if (cell_dict % has_key(id)) then + this % cells(i) = cell_dict % get_key(id) + else + call fatal_error("Could not find cell " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end do + + ! Generate mapping from cell indices to filter bins. + do i = 1, this % n_bins + call this % map % add_key(this % cells(i), i) + end do + end subroutine initialize_cell + + function text_label_cell(this, bin) result(label) + class(CellFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Cell " // to_str(cells(this % cells(bin)) % id) + end function text_label_cell + +end module tally_filter_cell diff --git a/src/tallies/tally_filter_cellborn.F90 b/src/tallies/tally_filter_cellborn.F90 new file mode 100644 index 000000000..35031ab00 --- /dev/null +++ b/src/tallies/tally_filter_cellborn.F90 @@ -0,0 +1,93 @@ +module tally_filter_cellborn + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error + use hdf5_interface + use geometry_header + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! CELLBORNFILTER specifies which cell the particle was born in. +!=============================================================================== + + type, public, extends(TallyFilter) :: CellbornFilter + integer, allocatable :: cells(:) + type(DictIntInt) :: map + contains + procedure :: get_all_bins => get_all_bins_cellborn + procedure :: to_statepoint => to_statepoint_cellborn + procedure :: text_label => text_label_cellborn + procedure :: initialize => initialize_cellborn + end type CellbornFilter + +contains + + subroutine get_all_bins_cellborn(this, p, estimator, match) + class(CellbornFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + if (this % map % has_key(p % cell_born)) then + call match % bins % push_back(this % map % get_key(p % cell_born)) + call match % weights % push_back(ONE) + end if + + end subroutine get_all_bins_cellborn + + subroutine to_statepoint_cellborn(this, filter_group) + class(CellbornFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: cell_ids(:) + + call write_dataset(filter_group, "type", "cellborn") + call write_dataset(filter_group, "n_bins", this % n_bins) + allocate(cell_ids(size(this % cells))) + do i = 1, size(this % cells) + cell_ids(i) = cells(this % cells(i)) % id + end do + call write_dataset(filter_group, "bins", cell_ids) + end subroutine to_statepoint_cellborn + + subroutine initialize_cellborn(this) + class(CellbornFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % cells(i) + if (cell_dict % has_key(id)) then + this % cells(i) = cell_dict % get_key(id) + else + call fatal_error("Could not find cell " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end do + + ! Generate mapping from cell indices to filter bins. + do i = 1, this % n_bins + call this % map % add_key(this % cells(i), i) + end do + end subroutine initialize_cellborn + + function text_label_cellborn(this, bin) result(label) + class(CellbornFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Birth Cell " // to_str(cells(this % cells(bin)) % id) + end function text_label_cellborn + +end module tally_filter_cellborn diff --git a/src/tallies/tally_filter_cellfrom.F90 b/src/tallies/tally_filter_cellfrom.F90 new file mode 100644 index 000000000..211743f2f --- /dev/null +++ b/src/tallies/tally_filter_cellfrom.F90 @@ -0,0 +1,77 @@ +module tally_filter_cellfrom + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error + use hdf5_interface + use geometry_header + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + use tally_filter_cell + + implicit none + private + +!=============================================================================== +! CELLFROMFILTER specifies which geometric cells particles exit when crossing a +! surface. +!=============================================================================== + + type, public, extends(CellFilter) :: CellFromFilter + contains + procedure :: get_all_bins => get_all_bins_cell_from + procedure :: to_statepoint => to_statepoint_cell_from + procedure :: text_label => text_label_cell_from + end type CellFromFilter + +contains + + subroutine get_all_bins_cell_from(this, p, estimator, match) + class(CellFromFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + + ! Starting one coordinate level deeper, find the next bin. + do i = 1, p % last_n_coord + if (this % map % has_key(p % last_cell(i))) then + call match % bins % push_back(this % map % get_key(p % last_cell(i))) + call match % weights % push_back(ONE) + exit + end if + end do + + end subroutine get_all_bins_cell_from + + subroutine to_statepoint_cell_from(this, filter_group) + class(CellFromFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: cell_ids(:) + + call write_dataset(filter_group, "type", "cellfrom") + call write_dataset(filter_group, "n_bins", this % n_bins) + + allocate(cell_ids(size(this % cells))) + do i = 1, size(this % cells) + cell_ids(i) = cells(this % cells(i)) % id + end do + call write_dataset(filter_group, "bins", cell_ids) + end subroutine to_statepoint_cell_from + + function text_label_cell_from(this, bin) result(label) + class(CellFromFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Cell from " // to_str(cells(this % cells(bin)) % id) + end function text_label_cell_from + +end module tally_filter_cellfrom diff --git a/src/tallies/tally_filter_delayedgroup.F90 b/src/tallies/tally_filter_delayedgroup.F90 new file mode 100644 index 000000000..d0ecf823a --- /dev/null +++ b/src/tallies/tally_filter_delayedgroup.F90 @@ -0,0 +1,60 @@ +module tally_filter_delayedgroup + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! DELAYEDGROUPFILTER bins outgoing fission neutrons in their delayed groups. +! The get_all_bins functionality is not actually used. The bins are manually +! iterated over in the scoring subroutines. +!=============================================================================== + + type, public, extends(TallyFilter) :: DelayedGroupFilter + integer, allocatable :: groups(:) + contains + procedure :: get_all_bins => get_all_bins_dg + procedure :: to_statepoint => to_statepoint_dg + procedure :: text_label => text_label_dg + end type DelayedGroupFilter + +contains + + subroutine get_all_bins_dg(this, p, estimator, match) + class(DelayedGroupFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + call match % bins % push_back(1) + call match % weights % push_back(ONE) + end subroutine get_all_bins_dg + + subroutine to_statepoint_dg(this, filter_group) + class(DelayedGroupFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "delayedgroup") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % groups) + end subroutine to_statepoint_dg + + function text_label_dg(this, bin) result(label) + class(DelayedGroupFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Delayed Group " // to_str(this % groups(bin)) + end function text_label_dg + +end module tally_filter_delayedgroup diff --git a/src/tallies/tally_filter_energyfunc.F90 b/src/tallies/tally_filter_energyfunc.F90 new file mode 100644 index 000000000..6e7245f20 --- /dev/null +++ b/src/tallies/tally_filter_energyfunc.F90 @@ -0,0 +1,91 @@ +module tally_filter_energyfunc + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use algorithm, only: binary_search + use constants + use error, only: fatal_error + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! EnergyFunctionFilter multiplies tally scores by an arbitrary function of +! incident energy described by a piecewise linear-linear interpolation. +!=============================================================================== + + type, public, extends(TallyFilter) :: EnergyFunctionFilter + real(8), allocatable :: energy(:) + real(8), allocatable :: y(:) + + contains + procedure :: get_all_bins => get_all_bins_energyfunction + procedure :: to_statepoint => to_statepoint_energyfunction + procedure :: text_label => text_label_energyfunction + end type EnergyFunctionFilter + +contains + + subroutine get_all_bins_energyfunction(this, p, estimator, match) + class(EnergyFunctionFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: n, indx + real(8) :: E, f, weight + + select type(this) + type is (EnergyFunctionFilter) + n = size(this % energy) + + ! Get pre-collision energy of particle + E = p % last_E + + ! Search to find incoming energy bin. + indx = binary_search(this % energy, n, E) + + ! Compute an interpolation factor between nearest bins. + f = (E - this % energy(indx)) & + / (this % energy(indx+1) - this % energy(indx)) + + ! Interpolate on the lin-lin grid. + call match % bins % push_back(1) + weight = (ONE - f) * this % y(indx) + f * this % y(indx+1) + call match % weights % push_back(weight) + end select + end subroutine get_all_bins_energyfunction + + subroutine to_statepoint_energyfunction(this, filter_group) + class(EnergyFunctionFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + select type(this) + type is (EnergyFunctionFilter) + call write_dataset(filter_group, "type", "energyfunction") + call write_dataset(filter_group, "energy", this % energy) + call write_dataset(filter_group, "y", this % y) + end select + end subroutine to_statepoint_energyfunction + + function text_label_energyfunction(this, bin) result(label) + class(EnergyFunctionFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + select type(this) + type is (EnergyFunctionFilter) + write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") & + "Energy Function f([", this % energy(1), ", ..., ", & + this % energy(size(this % energy)), "]) = [", this % y(1), & + ", ..., ", this % y(size(this % y)), "]" + end select + end function text_label_energyfunction + +end module tally_filter_energyfunc diff --git a/src/tallies/tally_filter_mu.F90 b/src/tallies/tally_filter_mu.F90 new file mode 100644 index 000000000..624dbb940 --- /dev/null +++ b/src/tallies/tally_filter_mu.F90 @@ -0,0 +1,74 @@ +module tally_filter_mu + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use algorithm, only: binary_search + use constants, only: ONE, MAX_LINE_LEN, NO_BIN_FOUND + use error, only: fatal_error + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! MUFILTER bins the incoming-outgoing direction cosine. This is only used for +! scatter reactions. +!=============================================================================== + + type, public, extends(TallyFilter) :: MuFilter + real(8), allocatable :: bins(:) + contains + procedure :: get_all_bins => get_all_bins_mu + procedure :: to_statepoint => to_statepoint_mu + procedure :: text_label => text_label_mu + end type MuFilter + +contains + + subroutine get_all_bins_mu(this, p, estimator, match) + class(MuFilter), 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 + + ! Search to find incoming energy bin. + bin = binary_search(this % bins, n + 1, p % mu) + if (bin /= NO_BIN_FOUND) then + call match % bins % push_back(bin) + call match % weights % push_back(ONE) + end if + end subroutine get_all_bins_mu + + subroutine to_statepoint_mu(this, filter_group) + class(MuFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "mu") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % bins) + end subroutine to_statepoint_mu + + function text_label_mu(this, bin) result(label) + class(MuFilter), 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 = "Change-in-Angle [" // trim(to_str(E0)) // ", " & + // trim(to_str(E1)) // ")" + end function text_label_mu + +end module tally_filter_mu diff --git a/src/tallies/tally_filter_polar.F90 b/src/tallies/tally_filter_polar.F90 new file mode 100644 index 000000000..9e3d9946b --- /dev/null +++ b/src/tallies/tally_filter_polar.F90 @@ -0,0 +1,82 @@ +module tally_filter_polar + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use algorithm, only: binary_search + use constants + use error, only: fatal_error + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! POLARFILTER bins the incident neutron polar angle (relative to the global +! z-axis). +!=============================================================================== + + type, public, extends(TallyFilter) :: PolarFilter + real(8), allocatable :: bins(:) + contains + procedure :: get_all_bins => get_all_bins_polar + procedure :: to_statepoint => to_statepoint_polar + procedure :: text_label => text_label_polar + end type PolarFilter + +contains + + subroutine get_all_bins_polar(this, p, estimator, match) + class(PolarFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: n + integer :: bin + real(8) :: theta + + n = this % n_bins + + ! Make sure the correct direction vector is used. + if (estimator == ESTIMATOR_TRACKLENGTH) then + theta = acos(p % coord(1) % uvw(3)) + else + theta = acos(p % last_uvw(3)) + end if + + ! Search to find polar angle bin. + bin = binary_search(this % bins, n + 1, theta) + if (bin /= NO_BIN_FOUND) then + call match % bins % push_back(bin) + call match % weights % push_back(ONE) + end if + end subroutine get_all_bins_polar + + subroutine to_statepoint_polar(this, filter_group) + class(PolarFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "polar") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % bins) + end subroutine to_statepoint_polar + + function text_label_polar(this, bin) result(label) + class(PolarFilter), 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 = "Polar Angle [" // trim(to_str(E0)) // ", " // trim(to_str(E1)) & + // ")" + end function text_label_polar + +end module tally_filter_polar diff --git a/src/tallies/tally_filter_surface.F90 b/src/tallies/tally_filter_surface.F90 new file mode 100644 index 000000000..aa8e390bf --- /dev/null +++ b/src/tallies/tally_filter_surface.F90 @@ -0,0 +1,92 @@ +module tally_filter_surface + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error + use hdf5_interface + use surface_header + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! SURFACEFILTER specifies which surface particles are crossing +!=============================================================================== + + type, public, extends(TallyFilter) :: SurfaceFilter + integer, allocatable :: surfaces(:) + + ! True if this filter is used for surface currents + logical :: current = .false. + contains + procedure :: get_all_bins => get_all_bins_surface + procedure :: to_statepoint => to_statepoint_surface + procedure :: text_label => text_label_surface + procedure :: initialize => initialize_surface + end type SurfaceFilter + +contains + + subroutine get_all_bins_surface(this, p, estimator, match) + class(SurfaceFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + + do i = 1, this % n_bins + if (abs(p % surface) == this % surfaces(i)) then + call match % bins % push_back(i) + if (p % surface < 0) then + call match % weights % push_back(-ONE) + else + call match % weights % push_back(ONE) + end if + exit + end if + end do + + end subroutine get_all_bins_surface + + subroutine to_statepoint_surface(this, filter_group) + class(SurfaceFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "surface") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", this % surfaces) + end subroutine to_statepoint_surface + + subroutine initialize_surface(this) + class(SurfaceFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % surfaces(i) + if (surface_dict % has_key(id)) then + this % surfaces(i) = surface_dict % get_key(id) + else + call fatal_error("Could not find surface " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end do + end subroutine initialize_surface + + function text_label_surface(this, bin) result(label) + class(SurfaceFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Surface " // to_str(surfaces(this % surfaces(bin)) % obj % id) + end function text_label_surface + +end module tally_filter_surface diff --git a/src/tallies/tally_filter_universe.F90 b/src/tallies/tally_filter_universe.F90 new file mode 100644 index 000000000..c48e4e917 --- /dev/null +++ b/src/tallies/tally_filter_universe.F90 @@ -0,0 +1,100 @@ +module tally_filter_universe + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error + use hdf5_interface + use geometry_header + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + +!=============================================================================== +! UNIVERSEFILTER specifies which geometric universes tally events reside in. +!=============================================================================== + + type, public, extends(TallyFilter) :: UniverseFilter + integer, allocatable :: universes(:) + type(DictIntInt) :: map + contains + procedure :: get_all_bins => get_all_bins_universe + procedure :: to_statepoint => to_statepoint_universe + procedure :: text_label => text_label_universe + procedure :: initialize => initialize_universe + end type UniverseFilter + +contains + + subroutine get_all_bins_universe(this, p, estimator, match) + class(UniverseFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + + ! Iterate over coordinate levels to see which universes match + do i = 1, p % n_coord + if (this % map % has_key(p % coord(i) % universe)) then + call match % bins % push_back(this % map % get_key(p % coord(i) & + % universe)) + call match % weights % push_back(ONE) + end if + end do + + end subroutine get_all_bins_universe + + subroutine to_statepoint_universe(this, filter_group) + class(UniverseFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: universe_ids(:) + + call write_dataset(filter_group, "type", "universe") + call write_dataset(filter_group, "n_bins", this % n_bins) + + allocate(universe_ids(size(this % universes))) + do i = 1, size(this % universes) + universe_ids(i) = universes(this % universes(i)) % id + end do + call write_dataset(filter_group, "bins", universe_ids) + end subroutine to_statepoint_universe + + subroutine initialize_universe(this) + class(UniverseFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % universes(i) + if (universe_dict % has_key(id)) then + this % universes(i) = universe_dict % get_key(id) + else + call fatal_error("Could not find universe " // trim(to_str(id)) & + &// " specified on a tally filter.") + end if + end do + + ! Generate mapping from universe indices to filter bins. + do i = 1, this % n_bins + call this % map % add_key(this % universes(i), i) + end do + end subroutine initialize_universe + + function text_label_universe(this, bin) result(label) + class(UniverseFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Universe " // to_str(universes(this % universes(bin)) % id) + end function text_label_universe + +end module tally_filter_universe