From 4d6972f2720064e202075204f102bcf35410cfe9 Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 00:17:22 -0400 Subject: [PATCH] only one score type for meshfilters and surface filters, updated docs, optimized last_cell data stored --- docs/source/io_formats/tallies.rst | 1 + docs/source/usersguide/tallies.rst | 28 ++--- src/geometry.F90 | 5 +- src/input_xml.F90 | 187 ++++++++++++++++------------- src/output.F90 | 2 +- src/particle_header.F90 | 26 +--- src/tally_filter.F90 | 6 +- src/tracking.F90 | 8 +- 8 files changed, 128 insertions(+), 135 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 62075be933..35a205a353 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -162,6 +162,7 @@ should be set to: surface IDs should be given. By default, net currents are tallied, and to tally a partial current from one cell to another, this should be used in combination with a cell or cell_from filter that defines the other cell. + This filter should not be used in combination with a meshfilter. :cellfrom: This filter allows the tally to be scored when crossing a surface and the diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 35b2ff9ac8..546dcabac2 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -261,21 +261,19 @@ The following tables show all valid scores: +----------------------+---------------------------------------------------+ |Score | Description | +======================+===================================================+ - |current |Partial currents on the boundaries of each cell in | - | |a mesh. Units are particles per source | - | |particle. Note that this score can only be used if | - | |a mesh filter has been specified. Furthermore, it | - | |may not be used in conjunction with any other | - | |score. Only energy and mesh filters may be used. | - +----------------------+---------------------------------------------------+ - |total-current |Total currents on any surface previously defined | - | |in the geometry. Units are particles per source | - | |particle. Note that this score cannot be used with | - | |a current score or a mesh filter. | - | |It may be used along with any other filter. | - | |Surfaces can be defined with cell_from and cell | - | |(to) filters or surface filters, partial currents | - | |can be obtained with cell_from and cell filters. | + |current |Used in combination with a mesh filter: | + | |Partial currents on the boundaries of each cell in | + | |a mesh. It may not be used in conjunction with any | + | |other score. Only energy and mesh filters may be | + | |used. | + | |Used in combination with a surface filter: | + | |Net currents on any surface previously defined in | + | |the geometry. It may be used along with any other | + | |filter, except mesh filters. | + | |Surfaces can alternatively be defined with cell | + | |from and cell filters and partial currents can be | + | |obtained by using cell_from and cell filters. | + | |Units are particles per source particle. | +----------------------+---------------------------------------------------+ |events |Number of scoring events. Units are events per | | |source particle. | diff --git a/src/geometry.F90 b/src/geometry.F90 index f80bd7eb70..c451b61278 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -384,9 +384,8 @@ contains ! the geometry, is reflected, or crosses into a new lattice or cell !=============================================================================== - subroutine cross_surface(p, last_cell) + subroutine cross_surface(p) type(Particle), intent(inout) :: p - integer, intent(in) :: last_cell ! last cell particle was in real(8) :: u ! x-component of direction real(8) :: v ! y-component of direction @@ -469,7 +468,7 @@ contains p%coord(1)%uvw(:) = [u, v, w] / norm ! Reassign particle's cell and surface - p % coord(1) % cell = last_cell + p % coord(1) % cell = p % last_cell(p % last_n_coord) p % surface = -p % surface ! If a reflective surface is coincident with a lattice or universe diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2bd6058b1f..e0fa63e0db 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3809,99 +3809,114 @@ contains case ('fission-q-recoverable') t % score_bins(j) = SCORE_FISS_Q_RECOV case ('current') - t % score_bins(j) = SCORE_CURRENT - t % type = TALLY_MESH_CURRENT - ! Check to make sure that current is the only desired response - ! for this tally - if (n_words > 1) then - call fatal_error("Cannot tally other scores in the & - &same tally as surface currents") - end if + ! Check which type of current is desired: mesh currents or + ! surface currents + if (t % find_filter(FILTER_SURFACE) > 0 .or. & + &t % find_filter(FILTER_CELL) > 0 .or. & + &t % find_filter(FILTER_CELLFROM) > 0) then - ! Get index of mesh filter - i_filter_mesh = t % filter(t % find_filter(FILTER_MESH)) - - ! Check to make sure mesh filter was specified - if (i_filter_mesh == 0) then - call fatal_error("Cannot tally surface current without a mesh & - &filter.") - end if - - ! Get pointer to mesh - select type(filt => filters(i_filter_mesh) % obj) - type is (MeshFilter) - i_mesh = filt % mesh - m => meshes(i_mesh) - end select - - ! Copy filter indices to temporary array - allocate(temp_filter(size(t % filter) + 1)) - temp_filter(1:size(t % filter)) = t % filter - - ! Move allocation back -- temp_filter becomes deallocated during - ! this call - call move_alloc(FROM=temp_filter, TO=t % filter) - n_filter = size(t % filter) - - ! Extend the filters array so we can add a surface filter and mesh - ! filter - call add_filters(2) - - ! Increment number of user filters - n_user_filters = n_user_filters + 2 - - ! Get index of the new mesh filter - i_filt = n_user_filters - 1 - - ! Duplicate the mesh filter since other tallies might use this - ! filter and we need to change the dimension - allocate(MeshFilter :: filters(i_filt) % obj) - select type(filt => filters(i_filt) % obj) - type is (MeshFilter) - filt % id = i_filt - filt % mesh = i_mesh - - ! We need to increase the dimension by one since we also need - ! currents coming into and out of the boundary mesh cells. - filt % n_bins = product(m % dimension + 1) - - ! Add filter to dictionary - call filter_dict % add_key(filt % id, i_filt) - end select - t % filter(t % find_filter(FILTER_MESH)) = i_filt - - ! Get index of the new surface filter - i_filt = n_user_filters - - ! Add surface filter - allocate(SurfaceFilter :: filters(i_filt) % obj) - select type (filt => filters(i_filt) % obj) - type is (SurfaceFilter) - filt % id = i_filt - filt % n_bins = 4 * m % n_dimension - allocate(filt % surfaces(4 * m % n_dimension)) - if (m % n_dimension == 1) then - filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT /) - elseif (m % n_dimension == 2) then - filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & - OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT /) - elseif (m % n_dimension == 3) then - filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & - OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT, OUT_BOTTOM, & - IN_BOTTOM, OUT_TOP, IN_TOP /) + ! Check to make sure that mesh currents are not desired as well + if (t % find_filter(FILTER_MESH) > 0) then + call fatal_error("Cannot tally other mesh currents & + &in the same tally as surface currents") end if - filt % current = .true. + + t % type = TALLY_SURFACE + t % score_bins(j) = SCORE_CURRENT + + else if (t % find_filter(FILTER_MESH) > 0) then + t % score_bins(j) = SCORE_CURRENT + t % type = TALLY_MESH_CURRENT + + ! Check to make sure that current is the only desired response + ! for this tally + if (n_words > 1) then + call fatal_error("Cannot tally other scores in the & + &same tally as surface currents") + end if + + ! Get index of mesh filter + i_filter_mesh = t % filter(t % find_filter(FILTER_MESH)) + + ! Check to make sure mesh filter was specified + if (i_filter_mesh == 0) then + call fatal_error("Cannot tally surface current without a mesh & + &filter.") + end if + + ! Get pointer to mesh + select type(filt => filters(i_filter_mesh) % obj) + type is (MeshFilter) + i_mesh = filt % mesh + m => meshes(i_mesh) + end select + + ! Copy filter indices to temporary array + allocate(temp_filter(size(t % filter) + 1)) + temp_filter(1:size(t % filter)) = t % filter + + ! Move allocation back -- temp_filter becomes deallocated during + ! this call + call move_alloc(FROM=temp_filter, TO=t % filter) + n_filter = size(t % filter) + + ! Extend the filters array so we can add a surface filter and + ! mesh filter + call add_filters(2) + + ! Increment number of user filters + n_user_filters = n_user_filters + 2 + + ! Get index of the new mesh filter + i_filt = n_user_filters - 1 + + ! Duplicate the mesh filter since other tallies might use this + ! filter and we need to change the dimension + allocate(MeshFilter :: filters(i_filt) % obj) + select type(filt => filters(i_filt) % obj) + type is (MeshFilter) + filt % id = i_filt + filt % mesh = i_mesh + + ! We need to increase the dimension by one since we also need + ! currents coming into and out of the boundary mesh cells. + filt % n_bins = product(m % dimension + 1) ! Add filter to dictionary call filter_dict % add_key(filt % id, i_filt) - end select - t % find_filter(FILTER_SURFACE) = n_filter - t % filter(n_filter) = i_filt + end select + t % filter(t % find_filter(FILTER_MESH)) = i_filt + + ! Get index of the new surface filter + i_filt = n_user_filters + + ! Add surface filter + allocate(SurfaceFilter :: filters(i_filt) % obj) + select type (filt => filters(i_filt) % obj) + type is (SurfaceFilter) + filt % id = i_filt + filt % n_bins = 4 * m % n_dimension + allocate(filt % surfaces(4 * m % n_dimension)) + if (m % n_dimension == 1) then + filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT /) + elseif (m % n_dimension == 2) then + filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & + OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT /) + elseif (m % n_dimension == 3) then + filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & + OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT, OUT_BOTTOM, & + IN_BOTTOM, OUT_TOP, IN_TOP /) + end if + filt % current = .true. + + ! Add filter to dictionary + call filter_dict % add_key(filt % id, i_filt) + end select + t % find_filter(FILTER_SURFACE) = n_filter + t % filter(n_filter) = i_filt + end if - case ('total-current') - t % type = TALLY_SURFACE - t % score_bins(j) = SCORE_CURRENT case ('events') t % score_bins(j) = SCORE_EVENTS case ('elastic', '(n,elastic)') diff --git a/src/output.F90 b/src/output.F90 index c9eecd24cb..ade1a5c05e 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -733,7 +733,7 @@ contains score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" - score_names(abs(SCORE_CURRENT)) = "Partial current" + score_names(abs(SCORE_CURRENT)) = "Current" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 7fcdf56683..76876f021d 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -49,8 +49,8 @@ module particle_header 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 + integer :: last_n_coord ! number of current coordinates + integer :: last_cell(MAX_COORD) ! coordinates for all levels ! Energy Data real(8) :: E ! post-collision energy @@ -110,7 +110,6 @@ module particle_header procedure :: clear => clear_particle procedure :: initialize_from_source procedure :: create_secondary - procedure :: set_last_coord end type Particle contains @@ -152,7 +151,6 @@ 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 @@ -224,26 +222,6 @@ 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 - this % last_coord = this % coord - - ! 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 - !=============================================================================== ! 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 1e40c409fc..908b6f1a65 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -746,7 +746,7 @@ contains start = 1 else do i = 1, p % last_n_coord - if (p % last_coord(i) % cell == this % cells(current_bin)) then + if (p % last_cell(i) == this % cells(current_bin)) then start = i + 1 exit end if @@ -757,8 +757,8 @@ contains next_bin = NO_BIN_FOUND weight = ERROR_REAL 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) + if (this % map % has_key(p % last_cell(i))) then + next_bin = this % map % get_key(p % last_cell(i)) weight = ONE exit end if diff --git a/src/tracking.F90 b/src/tracking.F90 index 8e67843e89..042e63c27f 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -150,8 +150,10 @@ contains if (next_level > 0) p % n_coord = next_level ! Saving previous cell data - last_cell = p % coord(p % n_coord) % cell - call p % set_last_coord() + do j = 1, p % n_coord + p % last_cell(j) = p % coord(j) % cell + end do + p % last_n_coord = p % n_coord ! Update last_ data. This is needed to use the same filters in ! surface tallies as the ones implemented for regular tallies @@ -168,7 +170,7 @@ contains ! Particle crosses surface p % surface = surface_crossed - call cross_surface(p, last_cell) + call cross_surface(p) p % event = EVENT_SURFACE end if ! Score cell to cell partial currents