crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #893 from GiudGiud/partial_proj

Browse source 
Partial current tallies as cell_to_cell tallies
This commit is contained in:
Sterling Harper 2017-08-08 18:44:55 -04:00 • committed by GitHub
commit 1ada0f4ec6
20 changed files with 801 additions and 137 deletions
Download patch file Download diff file Expand all files Collapse all files

29
docs/source/io_formats/tallies.rst
View file

@ -123,9 +123,9 @@ to the scored values. The ``filter`` element has the following
attributes/sub-elements:
:type:
The type of the filter. Accepted options are "cell", "cellborn",
"material", "universe", "energy", "energyout", "mu", "polar",
"azimuthal", "mesh", "distribcell", "delayedgroup", and
The type of the filter. Accepted options are "cell", "cellfrom",
"cellborn", "surface", "material", "universe", "energy", "energyout", "mu",
"polar", "azimuthal", "mesh", "distribcell", "delayedgroup", and
"energyfunction".
:bins:
@ -154,14 +154,31 @@ For each filter type, the following table describes what the ``bins`` attribute
should be set to:
:cell:
A list of unique IDs for cells in which the tally should be accumulated.
A list of unique IDs for cells in which the tally should be
accumulated.
:surface:
This filter allows the tally to be scored when crossing a surface. A list of
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
particle came from a specified cell. A list of cell IDs should be
given.
To tally a partial current from a cell to another, this filter should be
used in combination with a cell filter, to define the other cell.
This filter should not be used in combination with a meshfilter.
:cellborn:
This filter allows the tally to be scored to only when particles were
originally born in a specified cell. A list of cell IDs should be given.
originally born in a specified cell. A list of cell IDs should be
given.
:material:
A list of unique IDs for matreials in which the tally should be accumulated.
A list of unique IDs for materials in which the tally should be accumulated.
:universe:
A list of unique IDs for universes in which the tally should be accumulated.

1
docs/source/pythonapi/base.rst
View file

@ -117,6 +117,7 @@ Constructing Tallies
openmc.UniverseFilter
openmc.MaterialFilter
openmc.CellFilter
openmc.CellFromFilter
openmc.CellbornFilter
openmc.SurfaceFilter
openmc.MeshFilter

19
docs/source/usersguide/tallies.rst
View file

@ -261,12 +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. |
|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 thereby resulting in tallying|
| |partial currents. |
| |Units are particles per source particle. |
+----------------------+---------------------------------------------------+
|events |Number of scoring events. Units are events per |
| |source particle. |

35
openmc/filter.py
View file

@ -17,7 +17,7 @@ from .mixin import IDManagerMixin
_FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface',
'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal',
'distribcell', 'delayedgroup', 'energyfunction']
'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom']
_CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in',
3: 'x-max out', 4: 'x-max in',
@ -576,7 +576,40 @@ class CellFilter(WithIDFilter):
@property
def bins(self):
return self._bins
@bins.setter
def bins(self, bins):
self._smart_set_bins(bins, openmc.Cell)
class CellFromFilter(WithIDFilter):
"""Bins tally on which Cell the neutron came from.
Parameters
----------
bins : openmc.Cell, Integral, or iterable thereof
The Cell(s) to tally. Either openmc.Cell objects or their
Integral ID numbers can be used.
filter_id : int
Unique identifier for the filter
Attributes
----------
bins : Integral or Iterable of Integral
openmc.Cell IDs.
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
"""
@property
def bins(self):
return self._bins
@bins.setter
def bins(self, bins):
self._smart_set_bins(bins, openmc.Cell)

6
src/cmfd_input.F90
View file

@ -241,7 +241,7 @@ contains
! There are 3 tally types:
! 1: Only an energy in filter-> flux,total,p1 scatter
! 2: Energy in and energy out filter-> nu-scatter,nu-fission
! 3: Surface current
! 3: Mesh current
!===============================================================================
subroutine create_cmfd_tally(root)
@ -425,7 +425,7 @@ contains
end select
end if
! Duplicate the mesh filter for the surface current tally since other
! Duplicate the mesh filter for the mesh current tally since other
! tallies use this filter and we need to change the dimension
i_filt = i_filt + 1
allocate(MeshFilter :: filters(i_filt) % obj)
@ -597,7 +597,7 @@ contains
! Set macro bins
t % score_bins(1) = SCORE_CURRENT
t % type = TALLY_SURFACE_CURRENT
t % type = TALLY_MESH_CURRENT
end if
end do

10
src/constants.F90
View file

@ -288,7 +288,8 @@ module constants
! Tally type
integer, parameter :: &
TALLY_VOLUME = 1, &
TALLY_SURFACE_CURRENT = 2
TALLY_MESH_CURRENT = 2, &
TALLY_SURFACE = 3
! Tally estimator types
integer, parameter :: &
@ -318,7 +319,7 @@ module constants
SCORE_FISSION = -10, & ! fission rate
SCORE_NU_FISSION = -11, & ! neutron production rate
SCORE_KAPPA_FISSION = -12, & ! fission energy production rate
SCORE_CURRENT = -13, & ! partial current
SCORE_CURRENT = -13, & ! current
SCORE_FLUX_YN = -14, & ! angular moment of flux
SCORE_TOTAL_YN = -15, & ! angular moment of total reaction rate
SCORE_SCATTER_YN = -16, & ! angular flux-weighted scattering moment (0:N)
@ -352,7 +353,7 @@ module constants
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: N_FILTER_TYPES = 14
integer, parameter :: N_FILTER_TYPES = 15
integer, parameter :: &
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
@ -367,7 +368,8 @@ module constants
FILTER_POLAR = 11, &
FILTER_AZIMUTHAL = 12, &
FILTER_DELAYEDGROUP = 13, &
FILTER_ENERGYFUNCTION = 14
FILTER_ENERGYFUNCTION = 14, &
FILTER_CELLFROM = 15
! Mesh types
integer, parameter :: &

5
src/geometry.F90
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

2
src/global.F90
View file

@ -168,6 +168,7 @@ module global
type(VectorInt) :: active_current_tallies
type(VectorInt) :: active_collision_tallies
type(VectorInt) :: active_tallies
type(VectorInt) :: active_surface_tallies
! Global tallies
! 1) collision estimate of k-eff
@ -529,6 +530,7 @@ contains
call active_tracklength_tallies % clear()
call active_current_tallies % clear()
call active_collision_tallies % clear()
call active_surface_tallies % clear()
call active_tallies % clear()
! Deallocate track_identifiers

207
src/input_xml.F90
View file

@ -3007,7 +3007,7 @@ contains
end if
n_words = node_word_count(node_filt, "bins")
case ("mesh", "universe", "material", "cell", "distribcell", &
"cellborn", "surface", "delayedgroup")
"cellborn", "cellfrom", "surface", "delayedgroup")
if (.not. check_for_node(node_filt, "bins")) then
call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
end if
@ -3039,6 +3039,17 @@ contains
call get_node_array(node_filt, "bins", filt % cells)
end select
case ('cellfrom')
! Allocate and declare the filter type
allocate(CellFromFilter :: f % obj)
select type (filt => f % obj)
type is (CellFromFilter)
! Allocate and store bins
filt % n_bins = n_words
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
case ('cellborn')
! Allocate and declare the filter type
allocate(CellbornFilter :: f % obj)
@ -3073,7 +3084,6 @@ contains
end select
case ('surface')
call fatal_error("Surface filter is not yet supported!")
! Allocate and declare the filter type
allocate(SurfaceFilter :: f % obj)
select type (filt => f % obj)
@ -3404,6 +3414,8 @@ contains
t % find_filter(FILTER_DISTRIBCELL) = j
type is (CellFilter)
t % find_filter(FILTER_CELL) = j
type is (CellFromFilter)
t % find_filter(FILTER_CELLFROM) = j
type is (CellbornFilter)
t % find_filter(FILTER_CELLBORN) = j
type is (MaterialFilter)
@ -3441,13 +3453,6 @@ contains
! Store the filter indices
call move_alloc(FROM=temp_filter, TO=t % filter)
! Check that both cell and surface weren't specified
if (t % find_filter(FILTER_CELL) > 0 .and. &
t % find_filter(FILTER_SURFACE) > 0) then
call fatal_error("Cannot specify both cell and surface filters for &
&tally " // trim(to_str(t % id)))
end if
! =======================================================================
! READ DATA FOR NUCLIDES
@ -3816,95 +3821,113 @@ contains
case ('fission-q-recoverable')
t % score_bins(j) = SCORE_FISS_Q_RECOV
case ('current')
t % score_bins(j) = SCORE_CURRENT
t % type = TALLY_SURFACE_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 ('events')
t % score_bins(j) = SCORE_EVENTS

3
src/output.F90
View file

@ -758,6 +758,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)) = "Current"
! Create filename for tally output
filename = trim(path_output) // "tallies.out"
@ -816,7 +817,7 @@ contains
end if
! Handle surface current tallies separately
if (t % type == TALLY_SURFACE_CURRENT) then
if (t % type == TALLY_MESH_CURRENT) then
call write_surface_current(t, unit_tally)
cycle
end if

6
src/particle_header.F90
View file

@ -15,6 +15,7 @@ module particle_header
!===============================================================================
type LocalCoord
! Indices in various arrays for this level
integer :: cell = NONE
integer :: universe = NONE
@ -48,6 +49,10 @@ module particle_header
integer :: cell_instance ! offset for distributed properties
type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels
! Particle coordinates before crossing a surface
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
real(8) :: last_E ! pre-collision energy
@ -147,6 +152,7 @@ contains
! Set up base level coordinates
this % coord(1) % universe = root_universe
this % n_coord = 1
this % last_n_coord = 1
end subroutine initialize_particle

14
src/relaxng/tallies.rnc
View file

@ -35,12 +35,14 @@ element tallies {
element filter {
(element id { xsd:int } | attribute id { xsd:int }) &
(
( (element type { ( "cell" | "cellborn" | "material" | "universe" |
"surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" |
"polar" | "azimuthal" | "delayedgroup" | "energyfunction") } |
attribute type { ( "cell" | "cellborn" | "material" | "universe" |
"surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" |
"polar" | "azimuthal" | "delayedgroup" | "energyfunction") }) &
( (element type { ( "cell" | "cellfrom" | "cellborn" | "material" |
"universe" | "surface" | "distribcell" | "mesh" | "energy" |
"energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" |
"energyfunction") } |
attribute type { ( "cell" | "cellfrom" | "cellborn" | "material" |
"universe" | "surface" | "distribcell" | "mesh" | "energy" |
"energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" |
"energyfunction") }) &
(element bins { list { xsd:double+ } } |
attribute bins { list { xsd:double+ } })
) |

2
src/relaxng/tallies.rng
View file

@ -168,6 +168,7 @@
<element name="type">
<choice>
<value>cell</value>
<value>cellfrom</value>
<value>cellborn</value>
<value>material</value>
<value>universe</value>
@ -186,6 +187,7 @@
<attribute name="type">
<choice>
<value>cell</value>
<value>cellfrom</value>
<value>cellborn</value>
<value>material</value>
<value>universe</value>

147
src/tally.F90
View file

@ -3121,6 +3121,142 @@ contains
end subroutine score_collision_tally
!===============================================================================
! score_surface_tally is called at every surface crossing and can be used to
! tally total or partial currents between two cells
!===============================================================================
subroutine score_surface_tally(p)
type(Particle), intent(in) :: p
integer :: i
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: q ! loop index for scoring bins
integer :: k ! working index for expand and score
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
integer :: j ! loop index for scoring bins
integer :: filter_index ! single index for single bin
integer :: matching_bin ! next valid filter bin
real(8) :: flux ! collision estimate of flux
real(8) :: filter_weight ! combined weight of all filters
real(8) :: score ! analog tally score
logical :: finished ! found all valid bin combinations
! No collision, so no weight change when survival biasing
flux = p % wgt
TALLY_LOOP: do i = 1, active_surface_tallies % size()
! Get index of tally and pointer to tally
i_tally = active_surface_tallies % data(i)
associate (t => tallies(i_tally))
! Find all valid bins in each filter if they have not already been found
! for a previous tally.
do j = 1, size(t % filter)
i_filt = t % filter(j)
if (.not. filter_matches(i_filt) % bins_present) then
call filter_matches(i_filt) % bins % clear()
call filter_matches(i_filt) % weights % clear()
matching_bin = NO_BIN_FOUND
do
call filters(i_filt) % obj % get_next_bin(p, t % estimator, &
matching_bin, matching_bin, filter_weight)
if (matching_bin == NO_BIN_FOUND) exit
call filter_matches(i_filt) % bins % push_back(matching_bin)
call filter_matches(i_filt) % weights % push_back(filter_weight)
end do
filter_matches(i_filt) % bins_present = .true.
end if
! If there are no valid bins for this filter, then there is nothing to
! score and we can move on to the next tally.
if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
! Set the index of the bin used in the first filter combination
filter_matches(i_filt) % i_bin = 1
end do
! ========================================================================
! Loop until we've covered all valid bins on each of the filters.
FILTER_LOOP: do
! Reset scoring index and weight
filter_index = 1
filter_weight = ONE
! Determine scoring index and weight for this filter combination
do j = 1, size(t % filter)
i_filt = t % filter(j)
i_bin = filter_matches(i_filt) % i_bin
filter_index = filter_index + (filter_matches(i_filt) % bins % &
data(i_bin) - 1) * t % stride(j)
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
data(i_bin)
end do
! Determine score
score = flux * filter_weight
! Currently only one score type
k = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
k = k + 1
! determine what type of score bin
score_bin = t % score_bins(q)
! determine scoring bin index, no offset from nuclide bins
score_index = q
! Expand score if necessary and add to tally results.
call expand_and_score(p, t, score_index, filter_index, score_bin, &
score, k)
end do SCORE_LOOP
! ======================================================================
! Filter logic
! Increment the filter bins, starting with the last filter to find the
! next valid bin combination
finished = .true.
do j = size(t % filter), 1, -1
i_filt = t % filter(j)
if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
bins % size()) then
filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
finished = .false.
exit
else
filter_matches(i_filt) % i_bin = 1
end if
end do
! Once we have finished all valid bins for each of the filters, exit
! the loop.
if (finished) exit FILTER_LOOP
end do FILTER_LOOP
end associate
! If the user has specified that we can assume all tallies are spatially
! separate, this implies that once a tally has been scored to, we needn't
! check the others. This cuts down on overhead when there are many
! tallies specified
if (assume_separate) exit TALLY_LOOP
end do TALLY_LOOP
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
end subroutine score_surface_tally
!===============================================================================
! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually
! determining which mesh surfaces were crossed
@ -4353,9 +4489,12 @@ contains
elseif (user_tallies(i) % estimator == ESTIMATOR_COLLISION) then
call active_collision_tallies % push_back(i_user_tallies + i)
end if
elseif (user_tallies(i) % type == TALLY_SURFACE_CURRENT) then
elseif (user_tallies(i) % type == TALLY_MESH_CURRENT) then
call active_current_tallies % push_back(i_user_tallies + i)
elseif (user_tallies(i) % type == TALLY_SURFACE) then
call active_surface_tallies % push_back(i_user_tallies + i)
end if
end do
call active_tallies % shrink_to_fit()
@ -4363,6 +4502,7 @@ contains
call active_tracklength_tallies % shrink_to_fit()
call active_collision_tallies % shrink_to_fit()
call active_current_tallies % shrink_to_fit()
call active_surface_tallies % shrink_to_fit()
end subroutine setup_active_usertallies
@ -4386,6 +4526,9 @@ contains
else if (active_current_tallies % size() > 0) then
call fatal_error("Active current tallies should not exist before CMFD &
&tallies!")
else if (active_surface_tallies % size() > 0) then
call fatal_error("Active cell to cell tallies should not exist before &
&CMFD tallies!")
end if
do i = 1, n_cmfd_tallies
@ -4399,7 +4542,7 @@ contains
elseif (cmfd_tallies(i) % estimator == ESTIMATOR_TRACKLENGTH) then
call active_tracklength_tallies % push_back(i_cmfd_tallies + i)
end if
elseif (cmfd_tallies(i) % type == TALLY_SURFACE_CURRENT) then
elseif (cmfd_tallies(i) % type == TALLY_MESH_CURRENT) then
call active_current_tallies % push_back(i_cmfd_tallies + i)
end if
end do

82
src/tally_filter.F90
View file

@ -71,6 +71,17 @@ module tally_filter
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_next_bin => get_next_bin_cell_from
procedure :: to_statepoint => to_statepoint_cell_from
procedure :: text_label => text_label_cell_from
end type CellFromFilter
!===============================================================================
! DISTRIBCELLFILTER specifies which distributed geometric cells tally events
! reside in.
@ -98,8 +109,7 @@ module tally_filter
end type CellbornFilter
!===============================================================================
! SURFACEFILTER is currently not implemented for usual geometric surfaces, but
! it is used as a placeholder for mesh surfaces used in current tallies.
! SURFACEFILTER specifies which surface particles are crossing
!===============================================================================
type, extends(TallyFilter) :: SurfaceFilter
integer, allocatable :: surfaces(:)
@ -717,6 +727,70 @@ contains
label = "Cell " // to_str(cells(this % cells(bin)) % id)
end function text_label_cell
!===============================================================================
! CellFromFilter methods
!===============================================================================
subroutine get_next_bin_cell_from(this, p, estimator, current_bin, &
next_bin, weight)
class(CellFromFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
integer :: i, start
! 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_cell(i) == 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
do i = start, p % last_n_coord
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
end do
end subroutine get_next_bin_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
!===============================================================================
! DistribcellFilter methods
!===============================================================================
@ -887,9 +961,9 @@ contains
weight = ERROR_REAL
if (current_bin == NO_BIN_FOUND) then
do i = 1, this % n_bins
if (p % surface == this % surfaces(i)) then
if (abs(p % surface) == this % surfaces(i)) then
next_bin = i
weight = ONE
weight = sign(1, p % surface)
exit
end if
end do

21
src/tracking.F90
View file

@ -15,7 +15,7 @@ module tracking
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current, &
score_track_derivative, &
score_track_derivative, score_surface_tally, &
score_collision_derivative, zero_flux_derivs
use track_output, only: initialize_particle_track, write_particle_track, &
add_particle_track, finalize_particle_track
@ -36,7 +36,6 @@ contains
integer :: next_level ! next coordinate level to check
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 :: n_event ! number of collisions/crossings
real(8) :: d_boundary ! distance to nearest boundary
real(8) :: d_collision ! sampled distance to collision
@ -148,7 +147,18 @@ contains
! PARTICLE CROSSES SURFACE
if (next_level > 0) p % n_coord = next_level
last_cell = p % coord(p % n_coord) % cell
! Saving previous cell data
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
p % last_uvw = p % coord(p % n_coord) % uvw
p % last_E = p % E
p % coord(p % n_coord) % cell = NONE
if (any(lattice_translation /= 0)) then
! Particle crosses lattice boundary
@ -158,9 +168,12 @@ contains
else
! 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
if(active_surface_tallies % size() > 0) call score_surface_tally(p)
else
! ====================================================================
! PARTICLE HAS COLLISION

10
src/trigger.F90
View file

@ -167,8 +167,8 @@ contains
trigger % rel_err = ZERO
trigger % variance = ZERO
! Surface current tally triggers require special treatment
if (t % type == TALLY_SURFACE_CURRENT) then
! Mesh current tally triggers require special treatment
if (t % type == TALLY_MESH_CURRENT) then
call compute_tally_current(t, trigger)
else
@ -287,7 +287,7 @@ contains
!===============================================================================
! COMPUTE_TALLY_CURRENT computes the current for a surface current tally with
! COMPUTE_TALLY_CURRENT computes the current for a mesh current tally with
! precision trigger(s).
!===============================================================================
@ -307,8 +307,8 @@ contains
logical :: print_ebin ! should incoming energy bin be displayed?
real(8) :: rel_err = ZERO ! temporary relative error of result
real(8) :: std_dev = ZERO ! temporary standard deviration of result
type(TallyObject), pointer :: t ! surface current tally
type(TriggerObject) :: trigger ! surface current tally trigger
type(TallyObject), pointer :: t ! mesh current tally
type(TriggerObject) :: trigger ! mesh current tally trigger
type(RegularMesh), pointer :: m ! surface current mesh
! Get pointer to mesh

103
tests/test_surface_tally/inputs_true.dat Normal file
View file

@ -0,0 +1,103 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="13" material="13" name="fuel" region="-1" universe="9" />
<cell id="14" material="14" name="water" region="1 -3 9 -10" universe="9" />
<cell fill="9" id="15" name="root cell" region="2 -3 9 -10" universe="0" />
<surface coeffs="0 0 1" id="1" name="Fuel OR" type="z-cylinder" />
<surface boundary="vacuum" coeffs="-2" id="2" name="left" type="x-plane" />
<surface boundary="reflective" coeffs="2" id="3" name="right" type="x-plane" />
<surface boundary="reflective" coeffs="-2" id="9" name="bottom" type="y-plane" />
<surface boundary="reflective" coeffs="2" id="10" name="top" type="y-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="13" name="UO2 fuel at 2.4% wt enrichment">
<density units="g/cc" value="10.0" />
<nuclide ao="1.0" name="U238" />
<nuclide ao="0.02" name="U235" />
<nuclide ao="2.0" name="O16" />
</material>
<material id="14" name="Borated water">
<density units="g/cm3" value="1" />
<nuclide ao="0.0001" name="B10" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="fission">
<parameters>-0.62992 -0.62992 -1 0.62992 0.62992 1</parameters>
</space>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<filter id="5" type="cellfrom">
<bins>13</bins>
</filter>
<filter id="6" type="cell">
<bins>14</bins>
</filter>
<filter id="1" type="energy">
<bins>0.0 4000000.0 20000000.0</bins>
</filter>
<filter id="2" type="polar">
<bins>0.0 0.785398163397 3.14159265359</bins>
</filter>
<filter id="3" type="azimuthal">
<bins>0.0 0.785398163397 3.14159265359</bins>
</filter>
<filter id="4" type="surface">
<bins>1</bins>
</filter>
<filter id="7" type="cellfrom">
<bins>14</bins>
</filter>
<filter id="8" type="cell">
<bins>13</bins>
</filter>
<filter id="10" type="surface">
<bins>2</bins>
</filter>
<filter id="11" type="surface">
<bins>3</bins>
</filter>
<tally id="1" name="fuel_to_water_1">
<filters>5 6 1 2 3</filters>
<scores>current</scores>
</tally>
<tally id="2" name="fuel_to_water_2">
<filters>5 4 1 2 3</filters>
<scores>current</scores>
</tally>
<tally id="3" name="water_to_fuel_1">
<filters>7 8 1 2 3</filters>
<scores>current</scores>
</tally>
<tally id="4" name="water_to_fuel_2">
<filters>7 4 1 2 3</filters>
<scores>current</scores>
</tally>
<tally id="5" name="net_cylinder">
<filters>4 1 2 3</filters>
<scores>current</scores>
</tally>
<tally id="6" name="leakage_left">
<filters>10 1 2 3</filters>
<scores>current</scores>
</tally>
<tally id="7" name="net_right">
<filters>11 1</filters>
<scores>current</scores>
</tally>
<tally id="8" name="net_right">
<filters>11 1</filters>
<scores>current</scores>
</tally>
</tallies>

53
tests/test_surface_tally/results_true.dat Normal file
View file

@ -0,0 +1,53 @@
mean,std. dev.
2.1400000e-02,1.1850926e-03
7.4700000e-02,2.5649345e-03
1.6090000e-01,4.4433471e-03
6.4810000e-01,8.6041979e-03
2.0000000e-03,3.9440532e-04
4.2000000e-03,8.2731158e-04
1.1200000e-02,1.0934146e-03
4.6400000e-02,2.2715633e-03
2.1400000e-02,1.1850926e-03
7.4700000e-02,2.5649345e-03
1.6090000e-01,4.4433471e-03
6.4810000e-01,8.6041979e-03
2.0000000e-03,3.9440532e-04
4.2000000e-03,8.2731158e-04
1.1200000e-02,1.0934146e-03
4.6400000e-02,2.2715633e-03
5.5000000e-03,5.6273143e-04
4.4300000e-02,1.8502252e-03
4.5600000e-02,1.2220202e-03
4.1780000e-01,6.4460151e-03
0.0000000e+00,0.0000000e+00
1.4000000e-03,3.3993463e-04
0.0000000e+00,0.0000000e+00
1.6600000e-02,1.0561986e-03
-5.5000000e-03,5.6273143e-04
-4.4300000e-02,1.8502252e-03
-4.5600000e-02,1.2220202e-03
-4.1780000e-01,6.4460151e-03
0.0000000e+00,0.0000000e+00
-1.4000000e-03,3.3993463e-04
0.0000000e+00,0.0000000e+00
-1.6600000e-02,1.0561986e-03
1.5900000e-02,1.1200198e-03
3.0400000e-02,3.2734623e-03
1.1530000e-01,3.8327536e-03
2.3030000e-01,7.2250336e-03
2.0000000e-03,3.9440532e-04
2.8000000e-03,7.8598841e-04
1.1200000e-02,1.0934146e-03
2.9800000e-02,2.5811281e-03
0.0000000e+00,0.0000000e+00
-3.0000000e-02,1.3743685e-03
0.0000000e+00,0.0000000e+00
-3.4810000e-01,5.9711343e-03
0.0000000e+00,0.0000000e+00
-3.1000000e-03,6.9041051e-04
0.0000000e+00,0.0000000e+00
-2.1900000e-02,2.4241837e-03
-9.0250000e-01,2.0941320e-02
-3.5200000e-02,1.2806248e-03
-9.0250000e-01,2.0941320e-02
-3.5200000e-02,1.2806248e-03

183
tests/test_surface_tally/test_surface_tally.py Normal file
View file

@ -0,0 +1,183 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
import numpy as np
import openmc
import pandas as pd
class SurfaceTallyTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Instantiate some Materials and register the appropriate Nuclides
uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
uo2.set_density('g/cc', 10.0)
uo2.add_nuclide('U238', 1.0)
uo2.add_nuclide('U235', 0.02)
uo2.add_nuclide('O16', 2.0)
borated_water = openmc.Material(name='Borated water')
borated_water.set_density('g/cm3', 1)
borated_water.add_nuclide('B10', 10e-5)
borated_water.add_nuclide('H1', 2.0)
borated_water.add_nuclide('O16', 1.0)
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([uo2, borated_water])
materials_file.export_to_xml()
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=1, \
name='Fuel OR')
left = openmc.XPlane(surface_id=2, x0=-2, name='left')
right = openmc.XPlane(surface_id=3, x0=2, name='right')
bottom = openmc.YPlane(y0=-2, name='bottom')
top = openmc.YPlane(y0=2, name='top')
left.boundary_type = 'vacuum'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
# Instantiate Cells
fuel = openmc.Cell(name='fuel')
water = openmc.Cell(name='water')
# Use surface half-spaces to define regions
fuel.region = -fuel_or
water.region = +fuel_or & -right & +bottom & -top
# Register Materials with Cells
fuel.fill = uo2
water.fill = borated_water
# Instantiate pin cell Universe
pin_cell = openmc.Universe(name='pin cell')
pin_cell.add_cells([fuel, water])
# Instantiate root Cell and Universe
root_cell = openmc.Cell(name='root cell')
root_cell.region = +left & -right & +bottom & -top
root_cell.fill = pin_cell
root_univ = openmc.Universe(universe_id=0, name='root universe')
root_univ.add_cell(root_cell)
# Instantiate a Geometry, register the root Universe
geometry = openmc.Geometry(root_univ)
geometry.export_to_xml()
# Instantiate a Settings object, set all runtime parameters
settings_file = openmc.Settings()
settings_file.batches = 10
settings_file.inactive = 0
settings_file.particles = 1000
#settings_file.output = {'tallies': True}
# Create an initial uniform spatial source distribution
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:],\
only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
# Tallies file
tallies_file = openmc.Tallies()
# Create partial current tallies from fuel to water
# Filters
two_groups = [0., 4e6, 20e6]
energy_filter = openmc.EnergyFilter(two_groups)
polar_filter = openmc.PolarFilter([0, np.pi / 4, np.pi])
azimuthal_filter = openmc.AzimuthalFilter([0, np.pi / 4, np.pi])
surface_filter = openmc.SurfaceFilter([1])
cell_from_filter = openmc.CellFromFilter(fuel)
cell_filter = openmc.CellFilter(water)
# Use Cell to cell filters for partial current
cell_to_cell_tally = openmc.Tally(name=str('fuel_to_water_1'))
cell_to_cell_tally.filters = [cell_from_filter, cell_filter, \
energy_filter, polar_filter, azimuthal_filter]
cell_to_cell_tally.scores = ['current']
tallies_file.append(cell_to_cell_tally)
# Use a Cell from + surface filters for partial current
cell_to_cell_tally = openmc.Tally(name=str('fuel_to_water_2'))
cell_to_cell_tally.filters = [cell_from_filter, surface_filter, \
energy_filter, polar_filter, azimuthal_filter]
cell_to_cell_tally.scores = ['current']
tallies_file.append(cell_to_cell_tally)
# Create partial current tallies from water to fuel
# Filters
cell_from_filter = openmc.CellFromFilter(water)
cell_filter = openmc.CellFilter(fuel)
# Cell to cell filters for partial current
cell_to_cell_tally = openmc.Tally(name=str('water_to_fuel_1'))
cell_to_cell_tally.filters = [cell_from_filter, cell_filter, \
energy_filter, polar_filter, azimuthal_filter]
cell_to_cell_tally.scores = ['current']
tallies_file.append(cell_to_cell_tally)
# Cell from + surface filters for partial current
cell_to_cell_tally = openmc.Tally(name=str('water_to_fuel_2'))
cell_to_cell_tally.filters = [cell_from_filter, surface_filter, \
energy_filter, polar_filter, azimuthal_filter]
cell_to_cell_tally.scores = ['current']
tallies_file.append(cell_to_cell_tally)
# Create a net current tally on inner surface using a surface filter
surface_filter = openmc.SurfaceFilter([1])
surf_tally1 = openmc.Tally(name='net_cylinder')
surf_tally1.filters = [surface_filter, energy_filter, polar_filter, \
azimuthal_filter]
surf_tally1.scores = ['current']
tallies_file.append(surf_tally1)
# Create a net current tally on left surface using a surface filter
# This surface has a vacuum boundary condition, so leakage is tallied
surface_filter = openmc.SurfaceFilter([2])
surf_tally2 = openmc.Tally(name='leakage_left')
surf_tally2.filters = [surface_filter, energy_filter, polar_filter, \
azimuthal_filter]
surf_tally2.scores = ['current']
tallies_file.append(surf_tally2)
# Create a net current tally on right surface using a surface filter
# This surface has a reflective boundary condition, but the zero
# net current is not picked up because particles are only tallied once
surface_filter = openmc.SurfaceFilter([3])
surf_tally3 = openmc.Tally(name='net_right')
surf_tally3.filters = [surface_filter, energy_filter]
surf_tally3.scores = ['current']
tallies_file.append(surf_tally3)
surface_filter = openmc.SurfaceFilter([3])
surf_tally3 = openmc.Tally(name='net_right')
surf_tally3.filters = [surface_filter, energy_filter]
surf_tally3.scores = ['current']
tallies_file.append(surf_tally3)
tallies_file.export_to_xml()
def _get_results(self):
"""Digest info in the statepoint and return as a string."""
# Read the statepoint file.
sp = openmc.StatePoint(self._sp_name)
# Extract the tally data as a Pandas DataFrame.
df = pd.DataFrame()
for t in sp.tallies.values():
df = df.append(t.get_pandas_dataframe(), ignore_index=True)
# Extract the relevant data as a CSV string.
cols = ('mean', 'std. dev.')
return df.to_csv(None, columns=cols, index=False, float_format='%.7e')
return outstr
if __name__ == '__main__':
harness = SurfaceTallyTestHarness('statepoint.10.h5')
harness.main()