From d952d67fe545ba1d26a1a85fdcfa6206cfbfc9ae Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 08:39:16 -0400 Subject: [PATCH] made cell_from filter more general --- openmc/filter.py | 3 ++- src/geometry.F90 | 3 +-- src/particle_header.F90 | 55 ++++++++++++++++++++++++++++++++++++----- src/tally_filter.F90 | 49 ++++++++++++++++++++---------------- src/tracking.F90 | 10 +++----- 5 files changed, 83 insertions(+), 37 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index bc6b0bd35e..9382c333b3 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -18,7 +18,8 @@ AUTO_FILTER_ID = 10000 _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell', 'delayedgroup', 'energyfunction', 'celltocell', 'cellfrom', 'cellto'] + 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', + 'cellto'] _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in', 3: 'x-max out', 4: 'x-max in', diff --git a/src/geometry.F90 b/src/geometry.F90 index d6c6bff3b0..f80bd7eb70 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -11,8 +11,7 @@ module geometry use surface_header use stl_vector, only: VectorInt use string, only: to_str - use tally, only: score_surface_current, & - &score_partial_current + use tally, only: score_surface_current implicit none diff --git a/src/particle_header.F90 b/src/particle_header.F90 index bf22a2c6af..4f9aa25590 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -32,6 +32,7 @@ module particle_header logical :: rotated = .false. contains procedure :: reset => reset_coord + procedure :: copy_coord end type LocalCoord !=============================================================================== @@ -48,6 +49,10 @@ module particle_header integer :: n_coord ! number of current coordinates type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels + ! Particle coordinates before crossing a surface + integer :: last_n_coord ! number of current coordinates + type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels + ! Energy Data real(8) :: E ! post-collision energy real(8) :: last_E ! pre-collision energy @@ -59,10 +64,6 @@ module particle_header real(8) :: mu ! angle of scatter logical :: alive ! is particle alive? - ! Crossing surface data - integer :: last_cell ! last cell the particle was in before crossing - real(8) :: normal_proj ! cos of angle to surface normal when crossing - ! Pre-collision physical data real(8) :: last_xyz_current(3) ! coordinates of the last collision or ! reflective/periodic surface crossing @@ -110,6 +111,7 @@ module particle_header procedure :: clear => clear_particle procedure :: initialize_from_source procedure :: create_secondary + procedure :: set_last_coord end type Particle contains @@ -132,7 +134,6 @@ contains ! clear attributes this % surface = NONE - this % last_cell = NONE this % cell_born = NONE this % material = NONE this % last_material = NONE @@ -140,7 +141,6 @@ contains this % wgt = ONE this % last_wgt = ONE this % absorb_wgt = ZERO - this % normal_proj = 1 this % n_bank = 0 this % wgt_bank = ZERO this % sqrtkT = ERROR_REAL @@ -153,6 +153,8 @@ contains ! Set up base level coordinates this % coord(1) % universe = root_universe this % n_coord = 1 + this % last_coord(1) % universe = root_universe + this % last_n_coord = 1 end subroutine initialize_particle @@ -223,6 +225,47 @@ contains end subroutine initialize_from_source +!=============================================================================== +! SET_LAST_COORD copies all coordinate levels from coord to last_coord. This is +! meant to have information about the last cell in the cell_from filter +!=============================================================================== + + subroutine set_last_coord(this) + class(Particle), intent(inout) :: this + integer :: i + + this % last_n_coord = this % n_coord + + ! copy all information at all coordinate levels + do i = 1, this % n_coord + call this % last_coord(i) % copy_coord(this % coord(i)) + end do + + ! reset the rest of former last_coord + do i = this % n_coord + 1, MAX_COORD + call this % last_coord(i) % reset() + end do + + end subroutine set_last_coord + +!=============================================================================== +! COPY_COORD copies one coordinate levels from coord to last_coord. +!=============================================================================== + + elemental subroutine copy_coord(this, coord) + class(LocalCoord), intent(inout) :: this + class(LocalCoord), intent(in) :: coord + + this % cell = coord % cell + this % universe = coord % universe + this % lattice = coord % lattice + this % lattice_x = coord % lattice_x + this % lattice_y = coord % lattice_y + this % lattice_z = coord % lattice_z + this % rotated = coord % rotated + + end subroutine copy_coord + !=============================================================================== ! CREATE_SECONDARY stores the current phase space attributes of the particle in ! the secondary bank and increments the number of sites in the secondary bank. diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index eb36e8c9f0..c4e72eba5b 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -72,24 +72,21 @@ module tally_filter end type CellFilter !=============================================================================== -! CELLFROMFILTER specifies which geometric cells particles exit. +! CELLFROMFILTER specifies which geometric cells particles exit when crossing a +! surface. !=============================================================================== - type, extends(TallyFilter) :: CellFromFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map + type, extends(CellFilter) :: CellFromFilter contains - procedure :: get_next_bin => get_next_bin_cell_from - procedure :: to_statepoint => to_statepoint_cell_from - procedure :: text_label => text_label_cell_from - procedure :: initialize => initialize_cell_from + procedure :: get_next_bin_cell => get_next_bin_cell_from + procedure :: to_statepoint_cell => to_statepoint_cell_from + procedure :: text_label_cell => text_label_cell_from end type CellFromFilter !=============================================================================== -! CELLTOFILTER specifies which geometric cells particles enter. +! CELLTOFILTER specifies which geometric cells particles enter when crossing a +! surface. !=============================================================================== type, extends(CellFilter) :: CellToFilter -! integer, allocatable :: cells(:) -! type(DictIntInt) :: map contains procedure :: to_statepoint_cell => to_statepoint_cell_to procedure :: text_label_cell => text_label_cell_to @@ -755,18 +752,28 @@ contains integer :: i, start - ! Particle can only have one cell_from - ! If current_bin is already a bin, then no need to look for another bin + ! Find the coordinate level of the last bin we found. + if (current_bin == NO_BIN_FOUND) then + start = 1 + else + do i = 1, p % last_n_coord + if (p % last_coord(i) % cell == this % cells(current_bin)) then + start = i + 1 + exit + end if + end do + end if + + ! Starting one coordinate level deeper, find the next bin. next_bin = NO_BIN_FOUND weight = ERROR_REAL - if (current_bin == NO_BIN_FOUND) then - if (this % map % has_key(p % last_cell)) then - next_bin = this % map % get_key(p % last_cell) + do i = start, p % last_n_coord + if (this % map % has_key(p % last_coord(i) % cell)) then + next_bin = this % map % get_key(p % last_coord(i) % cell) weight = ONE + exit end if - else - next_bin = NO_BIN_FOUND - end if + end do end subroutine get_next_bin_cell_from @@ -777,7 +784,7 @@ contains integer :: i integer, allocatable :: cell_ids(:) - call write_dataset(filter_group, "type", "cell from") + call write_dataset(filter_group, "type", "cellfrom") call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) @@ -862,7 +869,7 @@ contains integer :: i integer, allocatable :: cell_ids(:) - call write_dataset(filter_group, "type", "cell to") + call write_dataset(filter_group, "type", "cellto") call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) diff --git a/src/tracking.F90 b/src/tracking.F90 index d8d0f8a7a5..c05cdf54fb 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -37,7 +37,6 @@ contains integer :: surface_crossed ! surface which particle is on integer :: lattice_translation(3) ! in-lattice translation vector integer :: last_cell ! most recent cell particle was in - integer :: last_n_coord ! most recent number of levels for last_cell integer :: n_event ! number of collisions/crossings real(8) :: d_boundary ! distance to nearest boundary real(8) :: d_collision ! sampled distance to collision @@ -152,6 +151,7 @@ contains ! Saving previous cell data last_cell = p % coord(p % n_coord) % cell +! p % set_last_coord() p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then @@ -163,9 +163,6 @@ contains ! Particle crosses surface p % surface = surface_crossed - ! Saving last cell for tallying purposes - p % last_cell = last_cell - ! Update last_ information. This is needed to use the same filters as ! the ones implemented for regular tallies p % last_uvw = p % coord(p % n_coord) % uvw @@ -173,10 +170,9 @@ contains call cross_surface(p, last_cell) p % event = EVENT_SURFACE - - ! Score cell to cell partial currents - call score_partial_current(p) end if + ! Score cell to cell partial currents + if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) else ! ==================================================================== ! PARTICLE HAS COLLISION