diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 5a5df4814d..691ecc09cc 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -55,6 +55,9 @@ _dll.openmc_material_filter_set_bins.errcheck = _error_handler _dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32] _dll.openmc_mesh_filter_set_mesh.restype = c_int _dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler +_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32] +_dll.openmc_meshsurface_filter_set_mesh.restype = c_int +_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler class Filter(_FortranObjectWithID): @@ -188,6 +191,10 @@ class MeshFilter(Filter): filter_type = 'mesh' +class MeshSurfaceFilter(Filter): + filter_type = 'meshsurface' + + class MuFilter(Filter): filter_type = 'mu' @@ -216,6 +223,7 @@ _FILTER_TYPE_MAP = { 'energyfunction': EnergyFunctionFilter, 'material': MaterialFilter, 'mesh': MeshFilter, + 'meshsurface': MeshSurfaceFilter, 'mu': MuFilter, 'polar': PolarFilter, 'surface': SurfaceFilter, diff --git a/src/api.F90 b/src/api.F90 index 64324da3b1..df183c292f 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -276,8 +276,9 @@ contains ! Clear active tally lists call active_analog_tallies % clear() call active_tracklength_tallies % clear() - call active_current_tallies % clear() + call active_meshsurf_tallies % clear() call active_collision_tallies % clear() + call active_surface_tallies % clear() call active_tallies % clear() ! Reset timers diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 5ab8971e6c..204284c480 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -35,14 +35,14 @@ element tallies { element filter { (element id { xsd:int } | attribute id { xsd:int }) & ( - ( (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 type { ( "cell" | "cellfrom" | "cellborn" | "material" | + "universe" | "surface" | "distribcell" | "mesh" | "energy" | + "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | + "energyfunction" | "meshsurface") } | + attribute type { ( "cell" | "cellfrom" | "cellborn" | "material" | + "universe" | "surface" | "distribcell" | "mesh" | "energy" | + "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | + "energyfunction" | "meshsurface") }) & (element bins { list { xsd:double+ } } | attribute bins { list { xsd:double+ } }) ) | diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index bad2fdce91..15b6f5b249 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -182,6 +182,7 @@ azimuthal delayedgroup energyfunction + meshsurface @@ -201,6 +202,7 @@ azimuthal delayedgroup energyfunction + meshsurface diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 058c414f22..f4d12be7e1 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3064,9 +3064,9 @@ contains ! tally total or partial currents between two cells !=============================================================================== - subroutine score_surface_tally(p) - - type(Particle), intent(in) :: p + subroutine score_surface_tally(p, tally_vec) + type(Particle), intent(in) :: p + type(VectorInt), intent(in) :: tally_vec integer :: i integer :: i_tally @@ -3086,9 +3086,9 @@ contains ! No collision, so no weight change when survival biasing flux = p % wgt - TALLY_LOOP: do i = 1, active_surface_tallies % size() + TALLY_LOOP: do i = 1, tally_vec % size() ! Get index of tally and pointer to tally - i_tally = active_surface_tallies % data(i) + i_tally = tally_vec % data(i) associate (t => tallies(i_tally) % obj) ! Find all valid bins in each filter if they have not already been found @@ -3188,134 +3188,6 @@ contains end subroutine score_surface_tally -!=============================================================================== -! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually -! determining which mesh surfaces were crossed -!=============================================================================== - - subroutine score_surface_current(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 - 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_current_tallies % size() - ! Get index of tally and pointer to tally - i_tally = active_current_tallies % data(i) - associate (t => tallies(i_tally) % obj) - - ! 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() - call filters(i_filt) % obj % get_all_bins(p, t % estimator, & - filter_matches(i_filt)) - 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 - - 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_current - !=============================================================================== ! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative !=============================================================================== @@ -4219,7 +4091,7 @@ contains call active_collision_tallies % clear() call active_tracklength_tallies % clear() call active_surface_tallies % clear() - call active_current_tallies % clear() + call active_meshsurf_tallies % clear() do i = 1, n_tallies associate (t => tallies(i) % obj) @@ -4237,7 +4109,7 @@ contains call active_collision_tallies % push_back(i) end if elseif (t % type == TALLY_MESH_SURFACE) then - call active_current_tallies % push_back(i) + call active_meshsurf_tallies % push_back(i) elseif (t % type == TALLY_SURFACE) then call active_surface_tallies % push_back(i) end if diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 6266005c2c..0495ce14dc 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -144,7 +144,7 @@ module tally_header ! Active tally lists type(VectorInt), public :: active_analog_tallies type(VectorInt), public :: active_tracklength_tallies - type(VectorInt), public :: active_current_tallies + type(VectorInt), public :: active_meshsurf_tallies type(VectorInt), public :: active_collision_tallies type(VectorInt), public :: active_tallies type(VectorInt), public :: active_surface_tallies @@ -418,7 +418,7 @@ contains ! Deallocate tally node lists call active_analog_tallies % clear() call active_tracklength_tallies % clear() - call active_current_tallies % clear() + call active_meshsurf_tallies % clear() call active_collision_tallies % clear() call active_surface_tallies % clear() call active_tallies % clear() diff --git a/src/tracking.F90 b/src/tracking.F90 index f5839eb8e9..3fe6a065b7 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -19,9 +19,9 @@ module tracking use surface_header use tally_header use tally, only: score_analog_tally, score_tracklength_tally, & - score_collision_tally, score_surface_current, & - score_track_derivative, score_surface_tally, & - score_collision_derivative, zero_flux_derivs + score_collision_tally, score_surface_tally, & + score_track_derivative, zero_flux_derivs, & + score_collision_derivative use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -185,7 +185,8 @@ contains p % event = EVENT_SURFACE end if ! Score cell to cell partial currents - if(active_surface_tallies % size() > 0) call score_surface_tally(p) + if(active_surface_tallies % size() > 0) & + call score_surface_tally(p, active_surface_tallies) else ! ==================================================================== ! PARTICLE HAS COLLISION @@ -200,7 +201,8 @@ contains ! since the direction of the particle will change and we need to use the ! pre-collision direction to figure out what mesh surfaces were crossed - if (active_current_tallies % size() > 0) call score_surface_current(p) + if (active_meshsurf_tallies % size() > 0) & + call score_surface_tally(p, active_meshsurf_tallies) ! Clear surface component p % surface = NONE @@ -317,12 +319,12 @@ contains ! forward slightly so that if the mesh boundary is on the surface, it is ! still processed - if (active_current_tallies % size() > 0) then + if (active_meshsurf_tallies % size() > 0) then ! TODO: Find a better solution to score surface currents than ! physically moving the particle forward slightly p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) + call score_surface_tally(p, active_meshsurf_tallies) end if ! Score to global leakage tally @@ -350,10 +352,10 @@ contains ! particle to change -- artificially move the particle slightly back in ! case the surface crossing is coincident with a mesh boundary - if (active_current_tallies % size() > 0) then + if (active_meshsurf_tallies % size() > 0) then xyz = p % coord(1) % xyz p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) + call score_surface_tally(p, active_meshsurf_tallies) p % coord(1) % xyz = xyz end if @@ -407,10 +409,10 @@ contains ! Score surface currents since reflection causes the direction of the ! particle to change -- artificially move the particle slightly back in ! case the surface crossing is coincident with a mesh boundary - if (active_current_tallies % size() > 0) then + if (active_meshsurf_tallies % size() > 0) then xyz = p % coord(1) % xyz p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) + call score_surface_tally(p, active_meshsurf_tallies) p % coord(1) % xyz = xyz end if diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index ca063b7ce9..511b298482 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -327,18 +327,27 @@ 1 + + 1 + 2 + + 2 + 3 + + 3 + 1 total - 1 + 4 current @@ -346,7 +355,7 @@ total - 2 + 5 current @@ -354,7 +363,7 @@ total - 3 + 6 current diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat index ac3439da50..f331238b70 100644 --- a/tests/regression_tests/filter_mesh/results_true.dat +++ b/tests/regression_tests/filter_mesh/results_true.dat @@ -1 +1 @@ -804d161cb8eae506d3247a533d122f44a01d3cedd566b3c65c71a0b51326dd9b97f8bbf45af7304b500476f3f854d91b10ccad94122d15f23641b05b835fada6 \ No newline at end of file +46950c046648faaa5ff3cb7b4fdd03667ae3c6da96a7ed8121761291de452cf88481f53e967ed52407d77d90109ae24839967a22ae216531a0e1491f5051ea72 \ No newline at end of file diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py index e0c6dd0f2f..8a7f7a6fca 100644 --- a/tests/regression_tests/filter_mesh/test.py +++ b/tests/regression_tests/filter_mesh/test.py @@ -30,6 +30,9 @@ class FilterMeshTestHarness(HashedPyAPITestHarness): mesh_1d_filter = openmc.MeshFilter(mesh_1d) mesh_2d_filter = openmc.MeshFilter(mesh_2d) mesh_3d_filter = openmc.MeshFilter(mesh_3d) + meshsurf_1d_filter = openmc.MeshSurfaceFilter(mesh_1d) + meshsurf_2d_filter = openmc.MeshSurfaceFilter(mesh_2d) + meshsurf_3d_filter = openmc.MeshSurfaceFilter(mesh_3d) # Initialized the tallies tally = openmc.Tally(name='tally 1') @@ -38,7 +41,7 @@ class FilterMeshTestHarness(HashedPyAPITestHarness): self._model.tallies.append(tally) tally = openmc.Tally(name='tally 2') - tally.filters = [mesh_1d_filter] + tally.filters = [meshsurf_1d_filter] tally.scores = ['current'] self._model.tallies.append(tally) @@ -48,7 +51,7 @@ class FilterMeshTestHarness(HashedPyAPITestHarness): self._model.tallies.append(tally) tally = openmc.Tally(name='tally 4') - tally.filters = [mesh_2d_filter] + tally.filters = [meshsurf_2d_filter] tally.scores = ['current'] self._model.tallies.append(tally) @@ -58,7 +61,7 @@ class FilterMeshTestHarness(HashedPyAPITestHarness): self._model.tallies.append(tally) tally = openmc.Tally(name='tally 6') - tally.filters = [mesh_3d_filter] + tally.filters = [meshsurf_3d_filter] tally.scores = ['current'] self._model.tallies.append(tally) diff --git a/tests/regression_tests/score_current/results_true.dat b/tests/regression_tests/score_current/results_true.dat index 75d61b7188..dedc473de9 100644 --- a/tests/regression_tests/score_current/results_true.dat +++ b/tests/regression_tests/score_current/results_true.dat @@ -1 +1 @@ -4675d5101f4f829369c39cb33d654430836b934ab07c165777ba6e214bdf3a698b8082f4f9bb9e78f1f0e495b30ea02cf9b3d14622c59915d818d678a1e5b7b1 \ No newline at end of file +138b312cdaa822c9b62f757b3259522004b679c4ed289a92798b29a6442c26d12c53256635be273f13e3703816ff50a1b9f52d79770eade01e482384ac0d389f \ No newline at end of file diff --git a/tests/regression_tests/score_current/tallies.xml b/tests/regression_tests/score_current/tallies.xml index 3b8f60adba..a381e4c2e2 100644 --- a/tests/regression_tests/score_current/tallies.xml +++ b/tests/regression_tests/score_current/tallies.xml @@ -9,7 +9,7 @@ - mesh + meshsurface 1