diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 9c0a2515e1..35a205a353 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -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. diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index dde1c91199..e6cf1e1ccc 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -117,6 +117,7 @@ Constructing Tallies openmc.UniverseFilter openmc.MaterialFilter openmc.CellFilter + openmc.CellFromFilter openmc.CellbornFilter openmc.SurfaceFilter openmc.MeshFilter diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 3ee2e99d71..3a742b59e8 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -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. | diff --git a/openmc/filter.py b/openmc/filter.py index 90f749668d..2b329ef0e0 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -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) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index d57c21beb8..68cd94aa6b 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -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 diff --git a/src/constants.F90 b/src/constants.F90 index 16c02d01b1..73aca25d29 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 :: & diff --git a/src/geometry.F90 b/src/geometry.F90 index 4444562321..7f75a817f6 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/global.F90 b/src/global.F90 index 2f47a9dce0..8343b49ac4 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2fc6ef9188..efdc3e5230 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index 9835225e06..b239a4c79b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/particle_header.F90 b/src/particle_header.F90 index f28e549761..3f4680c2ff 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -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 diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 4674fd4c03..5ab8971e6c 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -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+ } }) ) | diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 1c68e57e83..bad2fdce91 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -168,6 +168,7 @@ cell + cellfrom cellborn material universe @@ -186,6 +187,7 @@ cell + cellfrom cellborn material universe diff --git a/src/tally.F90 b/src/tally.F90 index cf62412f03..1bdde1b506 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index dfb89d4330..908b6f1a65 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -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 diff --git a/src/tracking.F90 b/src/tracking.F90 index 7650f598ab..546ae38809 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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 diff --git a/src/trigger.F90 b/src/trigger.F90 index dc0a10b153..cb5f2c4c85 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -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 diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat new file mode 100644 index 0000000000..e66d44273a --- /dev/null +++ b/tests/test_surface_tally/inputs_true.dat @@ -0,0 +1,103 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 0 + + + -0.62992 -0.62992 -1 0.62992 0.62992 1 + + + + + + + 13 + + + 14 + + + 0.0 4000000.0 20000000.0 + + + 0.0 0.785398163397 3.14159265359 + + + 0.0 0.785398163397 3.14159265359 + + + 1 + + + 14 + + + 13 + + + 2 + + + 3 + + + 5 6 1 2 3 + current + + + 5 4 1 2 3 + current + + + 7 8 1 2 3 + current + + + 7 4 1 2 3 + current + + + 4 1 2 3 + current + + + 10 1 2 3 + current + + + 11 1 + current + + + 11 1 + current + + diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat new file mode 100644 index 0000000000..51e2d7dc0a --- /dev/null +++ b/tests/test_surface_tally/results_true.dat @@ -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 diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py new file mode 100644 index 0000000000..ae525a6838 --- /dev/null +++ b/tests/test_surface_tally/test_surface_tally.py @@ -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()