only one score type for meshfilters and surface filters, updated docs, optimized last_cell data stored

This commit is contained in:
guillaume 2017-07-14 00:17:22 -04:00
parent 9fd8bd3fcf
commit 4d6972f272
8 changed files with 128 additions and 135 deletions

View file

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

View file

@ -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. |

View file

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

View file

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

View file

@ -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"

View file

@ -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.

View file

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

View file

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