made cell_from filter more general

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',

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
 

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.

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

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