From 794ac9255ae077a1dc6e496f24df194d5d669099 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 15 Mar 2018 14:07:42 -0500 Subject: [PATCH] Add initial Fortran implementation of MeshSurfaceFilter --- CMakeLists.txt | 1 + src/constants.F90 | 5 +- src/input_xml.F90 | 11 +- src/tallies/tally_filter.F90 | 5 + src/tallies/tally_filter_mesh.F90 | 3 - src/tallies/tally_filter_meshsurface.F90 | 392 +++++++++++++++++++++++ src/tallies/tally_header.F90 | 2 + 7 files changed, 411 insertions(+), 8 deletions(-) create mode 100644 src/tallies/tally_filter_meshsurface.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index be62dcd57b..9517fe2ec3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -422,6 +422,7 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/tally_filter_energyfunc.F90 src/tallies/tally_filter_material.F90 src/tallies/tally_filter_mesh.F90 + src/tallies/tally_filter_meshsurface.F90 src/tallies/tally_filter_mu.F90 src/tallies/tally_filter_polar.F90 src/tallies/tally_filter_surface.F90 diff --git a/src/constants.F90 b/src/constants.F90 index 3334a4cc76..75a9471490 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -358,7 +358,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 15 + integer, parameter :: N_FILTER_TYPES = 16 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -374,7 +374,8 @@ module constants FILTER_AZIMUTHAL = 12, & FILTER_DELAYEDGROUP = 13, & FILTER_ENERGYFUNCTION = 14, & - FILTER_CELLFROM = 15 + FILTER_CELLFROM = 15, & + FILTER_MESHSURFACE = 16 ! Mesh types integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 918d471d3c..792450e517 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -10,7 +10,7 @@ module input_xml use distribution_multivariate use distribution_univariate use endf, only: reaction_name - use error, only: fatal_error, warning, write_message + use error, only: fatal_error, warning, write_message, openmc_err_msg use geometry, only: calc_offsets, maximum_levels, count_instance, & neighbor_lists use geometry_header @@ -2358,13 +2358,13 @@ contains call get_node_value(node_filt, "type", temp_str) temp_str = to_lower(temp_str) - ! Determine number of bins + ! Make sure bins have been set select case(temp_str) case ("energy", "energyout", "mu", "polar", "azimuthal") if (.not. check_for_node(node_filt, "bins")) then call fatal_error("Bins not set in filter " // trim(to_str(filter_id))) end if - case ("mesh", "universe", "material", "cell", "distribcell", & + case ("mesh", "meshsurface", "universe", "material", "cell", "distribcell", & "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))) @@ -2373,6 +2373,9 @@ contains ! Allocate according to the filter type err = openmc_filter_set_type(i_start + i - 1, to_c_string(temp_str)) + if (err /= 0) then + call fatal_error(to_f_string(openmc_err_msg)) + end if ! Read filter data from XML call f % obj % from_xml(node_filt) @@ -2923,6 +2926,8 @@ contains ! Set filters err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter) deallocate(temp_filter) + else + t % score_bins(j) = SCORE_CURRENT end if case ('events') diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 587f8dd902..3d6eaaef21 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -19,6 +19,7 @@ module tally_filter use tally_filter_energyfunc use tally_filter_material use tally_filter_mesh + use tally_filter_meshsurface use tally_filter_mu use tally_filter_polar use tally_filter_surface @@ -67,6 +68,8 @@ contains type_ = 'material' type is (MeshFilter) type_ = 'mesh' + type is (MeshSurfaceFilter) + type_ = 'meshsurface' type is (MuFilter) type_ = 'mu' type is (PolarFilter) @@ -136,6 +139,8 @@ contains allocate(MaterialFilter :: filters(index) % obj) case ('mesh') allocate(MeshFilter :: filters(index) % obj) + case ('meshsurface') + allocate(MeshSurfaceFilter :: filters(index) % obj) case ('mu') allocate(MuFilter :: filters(index) % obj) case ('polar') diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 2092e7e717..3c0e6250c2 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -84,7 +84,6 @@ contains integer :: ijk1(3) ! indices of ending coordinates integer :: search_iter ! loop count for intersection search integer :: bin - real(8) :: weight ! weight to be pushed back real(8) :: uvw(3) ! cosine of angle of particle real(8) :: xyz0(3) ! starting/intermediate coordinates real(8) :: xyz1(3) ! ending coordinates of particle @@ -96,8 +95,6 @@ contains logical :: end_in_mesh ! ending coordinates inside mesh? type(RegularMesh), pointer :: m - weight = ERROR_REAL - ! Get a pointer to the mesh. m => meshes(this % mesh) n = m % n_dimension diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 new file mode 100644 index 0000000000..e45261f440 --- /dev/null +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -0,0 +1,392 @@ +module tally_filter_meshsurface + + use, intrinsic :: ISO_C_BINDING + + use hdf5 + + use constants + use dict_header, only: EMPTY + use error + use mesh_header, only: RegularMesh, meshes, n_meshes, mesh_dict + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + use xml_interface + + implicit none + private + public :: openmc_meshsurface_filter_set_mesh + +!=============================================================================== +! MESHFILTER indexes the location of particle events to a regular mesh. For +! tracklength tallies, it will produce multiple valid bins and the bin weight +! will correspond to the fraction of the track length that lies in that bin. +!=============================================================================== + + type, public, extends(TallyFilter) :: MeshSurfaceFilter + integer :: mesh + contains + procedure :: from_xml + procedure :: get_all_bins + procedure :: to_statepoint + procedure :: text_label + end type MeshSurfaceFilter + +contains + + subroutine from_xml(this, node) + class(MeshSurfaceFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + + integer :: i_mesh + integer :: id + integer :: n + integer :: n_dim + integer :: val + + n = node_word_count(node, "bins") + + if (n /= 1) call fatal_error("Only one mesh can be & + &specified per mesh filter.") + + ! Determine id of mesh + call get_node_value(node, "bins", id) + + ! Get pointer to mesh + val = mesh_dict % get(id) + if (val /= EMPTY) then + i_mesh = val + else + call fatal_error("Could not find mesh " // trim(to_str(id)) & + // " specified on filter.") + end if + + ! Determine number of bins + n_dim = meshes(i_mesh) % n_dimension + this % n_bins = 4*n_dim*product(meshes(i_mesh) % dimension + 1) + + ! Store the index of the mesh + this % mesh = i_mesh + end subroutine from_xml + + subroutine get_all_bins(this, p, estimator, match) + class(MeshSurfaceFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: j ! loop indices + integer :: n_dim ! num dimensions of the mesh + integer :: d1 ! dimension index + integer :: d2 ! dimension index + integer :: d3 ! dimension index + integer :: ijk0(3) ! indices of starting coordinates + integer :: ijk1(3) ! indices of ending coordinates + integer :: n_cross ! number of surface crossings + integer :: i_mesh ! flattened mesh bin index + integer :: i_surf ! surface index (1--12) + integer :: i_bin ! actual index for filter + real(8) :: uvw(3) ! cosine of angle of particle + real(8) :: xyz0(3) ! starting/intermediate coordinates + real(8) :: xyz1(3) ! ending coordinates of particle + real(8) :: xyz_cross(3) ! coordinates of bounding surfaces + real(8) :: d(3) ! distance to each bounding surface + real(8) :: distance ! actual distance traveled + logical :: start_in_mesh ! particle's starting xyz in mesh? + logical :: end_in_mesh ! particle's ending xyz in mesh? + logical :: cross_surface ! whether the particle crosses a surface + + ! Copy starting and ending location of particle + xyz0 = p % last_xyz_current + xyz1 = p % coord(1) % xyz + + associate (m => meshes(this % mesh)) + n_dim = m % n_dimension + + ! Determine indices for starting and ending location + call m % get_indices(xyz0, ijk0, start_in_mesh) + call m % get_indices(xyz1, ijk1, end_in_mesh) + + ! Check to see if start or end is in mesh -- if not, check if track still + ! intersects with mesh + if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + if (.not. m % intersects(xyz0, xyz1)) return + end if + + ! Calculate number of surface crossings + n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim))) + if (n_cross == 0) return + + ! Copy particle's direction + uvw = p % coord(1) % uvw + + ! Bounding coordinates + do d1 = 1, n_dim + if (uvw(d1) > 0) then + xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) + else + xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) + end if + end do + + do j = 1, n_cross + ! Set the distances to infinity + d = INFINITY + + ! Calculate distance to each bounding surface. We need to treat + ! special case where the cosine of the angle is zero since this would + ! result in a divide-by-zero. + do d1 = 1, n_dim + if (uvw(d1) == 0) then + d(d1) = INFINITY + else + d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) + end if + end do + + ! Determine the closest bounding surface of the mesh cell by + ! calculating the minimum distance. Then use the minimum distance and + ! direction of the particle to determine which surface was crossed. + distance = minval(d) + + ! Loop over the dimensions + do d1 = 1, n_dim + + ! Get the other dimensions. + if (d1 == 1) then + d2 = mod(d1, 3) + 1 + d3 = mod(d1 + 1, 3) + 1 + else + d2 = mod(d1 + 1, 3) + 1 + d3 = mod(d1, 3) + 1 + end if + + ! Check whether distance is the shortest distance + if (distance == d(d1)) then + + ! Check whether particle is moving in positive d1 direction + if (uvw(d1) > 0) then + + ! Outward current on d1 max surface + if (all(ijk0(:n_dim) >= 1) .and. & + all(ijk0(:n_dim) <= m % dimension)) then + i_surf = d1 * 4 - 1 + i_mesh = m % get_bin_from_indices(ijk0) + i_bin = 4*n_dim*i_mesh + i_surf + 1 + + call match % bins % push_back(i_bin) + call match % weights % push_back(ONE) + end if + + ! Inward current on d1 min surface + cross_surface = .false. + select case(n_dim) + + case (1) + if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1)) then + cross_surface = .true. + end if + + case (2) + if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then + cross_surface = .true. + end if + + case (3) + if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then + cross_surface = .true. + end if + end select + + ! If the particle crossed the surface, tally the current + if (cross_surface) then + ijk0(d1) = ijk0(d1) + 1 + i_surf = d1 * 4 - 2 + i_mesh = m % get_bin_from_indices(ijk0) + i_bin = 4*n_dim*i_mesh + i_surf + 1 + + call match % bins % push_back(i_bin) + call match % weights % push_back(ONE) + + ijk0(d1) = ijk0(d1) - 1 + end if + + ijk0(d1) = ijk0(d1) + 1 + xyz_cross(d1) = xyz_cross(d1) + m % width(d1) + + ! The particle is moving in the negative d1 direction + else + + ! Outward current on d1 min surface + if (all(ijk0(:n_dim) >= 1) .and. & + all(ijk0(:n_dim) <= m % dimension)) then + i_surf = d1 * 4 - 3 + i_mesh = m % get_bin_from_indices(ijk0) + i_bin = 4*n_dim*i_mesh + i_surf + 1 + + call match % bins % push_back(i_bin) + call match % weights % push_back(ONE) + end if + + ! Inward current on d1 max surface + cross_surface = .false. + select case(n_dim) + + case (1) + if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1) then + cross_surface = .true. + end if + + case (2) + if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.& + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then + cross_surface = .true. + end if + + case (3) + if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.& + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then + cross_surface = .true. + end if + end select + + ! If the particle crossed the surface, tally the current + if (cross_surface) then + ijk0(d1) = ijk0(d1) - 1 + i_surf = d1 * 4 + i_mesh = m % get_bin_from_indices(ijk0) + i_bin = 4*n_dim*i_mesh + i_surf + 1 + + call match % bins % push_back(i_bin) + call match % weights % push_back(ONE) + + ijk0(d1) = ijk0(d1) + 1 + end if + + ijk0(d1) = ijk0(d1) - 1 + xyz_cross(d1) = xyz_cross(d1) - m % width(d1) + end if + end if + end do + + ! Calculate new coordinates + xyz0 = xyz0 + distance * uvw + end do + end associate + + end subroutine get_all_bins + + subroutine to_statepoint(this, filter_group) + class(MeshSurfaceFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "meshsurface") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", meshes(this % mesh) % id) + end subroutine to_statepoint + + function text_label(this, bin) result(label) + class(MeshSurfaceFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + integer :: i_mesh + integer :: i_surf + integer :: n_dim + integer, allocatable :: ijk(:) + + associate (m => meshes(this % mesh)) + n_dim = m % n_dimension + allocate(ijk(n_dim)) + + ! Get flattend mesh index and surface index + i_mesh = (bin - 1) / (4*n_dim) + i_surf = mod(bin - 1, 4*n_dim) + 1 + + ! Get mesh index part of label + call m % get_indices_from_bin(i_mesh, ijk) + if (m % n_dimension == 1) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" + elseif (m % n_dimension == 2) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" + elseif (m % n_dimension == 3) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + end if + + ! Get surface part of label + select case (i_surf) + case (OUT_LEFT) + label = trim(label) // " Outgoing, x-min" + case (IN_LEFT) + label = trim(label) // " Incoming, x-min" + case (OUT_RIGHT) + label = trim(label) // " Outgoing, x-max" + case (IN_RIGHT) + label = trim(label) // " Incoming, x-max" + case (OUT_BACK) + label = trim(label) // " Outgoing, y-min" + case (IN_BACK) + label = trim(label) // " Incoming, y-min" + case (OUT_FRONT) + label = trim(label) // " Outgoing, y-max" + case (IN_FRONT) + label = trim(label) // " Incoming, y-max" + case (OUT_BOTTOM) + label = trim(label) // " Outgoing, z-min" + case (IN_BOTTOM) + label = trim(label) // " Incoming, z-min" + case (OUT_TOP) + label = trim(label) // " Outgoing, z-max" + case (IN_TOP) + label = trim(label) // " Incoming, z-max" + end select + end associate + end function text_label + +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) bind(C) + ! Set the mesh for a mesh filter + integer(C_INT32_T), value, intent(in) :: index + integer(C_INT32_T), value, intent(in) :: index_mesh + integer(C_INT) :: err + + integer :: n_dim + + err = 0 + if (index >= 1 .and. index <= n_filters) then + if (allocated(filters(index) % obj)) then + select type (f => filters(index) % obj) + type is (MeshSurfaceFilter) + if (index_mesh >= 1 .and. index_mesh <= n_meshes) then + f % mesh = index_mesh + n_dim = meshes(index_mesh) % n_dimension + f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension + 1) + else + err = E_OUT_OF_BOUNDS + call set_errmsg("Index in 'meshes' array is out of bounds.") + end if + class default + err = E_INVALID_TYPE + call set_errmsg("Tried to set mesh on a non-mesh filter.") + end select + else + err = E_ALLOCATE + call set_errmsg("Filter type has not been set yet.") + end if + else + err = E_OUT_OF_BOUNDS + call set_errmsg("Index in filters array out of bounds.") + end if + end function openmc_meshsurface_filter_set_mesh + +end module tally_filter_meshsurface diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 9604d2382a..6266005c2c 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -336,6 +336,8 @@ contains j = FILTER_SURFACE type is (MeshFilter) j = FILTER_MESH + type is (MeshSurfaceFilter) + j = FILTER_MESHSURFACE type is (EnergyFilter) j = FILTER_ENERGYIN type is (EnergyoutFilter)