diff --git a/src/cell.cpp b/src/cell.cpp index c22813ea7..a3a5d9fde 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -692,18 +692,6 @@ extern "C" { double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT[i];} - bool cell_simple(Cell* c) {return c->simple;} - - void cell_distance(Cell* c, double xyz[3], double uvw[3], int32_t on_surface, - double* min_dist, int32_t* i_surf) - { - Position r {xyz}; - Direction u {uvw}; - std::pair out = c->distance(r, u, on_surface); - *min_dist = out.first; - *i_surf = out.second; - } - int32_t cell_offset(Cell* c, int map) {return c->offset[map];} void cell_to_hdf5(Cell* c, hid_t group) {c->to_hdf5(group);} diff --git a/src/geometry.F90 b/src/geometry.F90 index dcd44360b..b5dee81fa 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -43,6 +43,16 @@ module geometry integer(C_INT), intent(in) :: lattice_translation(3) end subroutine cross_lattice + subroutine distance_to_boundary(p, dist, surface_crossed, & + lattice_translation, next_level) bind(C) + import Particle, C_DOUBLE, C_INT + type(Particle), intent(inout) :: p + real(C_DOUBLE), intent(out) :: dist + integer(C_INT), intent(out) :: surface_crossed + integer(C_INT), intent(out) :: lattice_translation(3) + integer(C_INT), intent(out) :: next_level + end subroutine distance_to_boundary + subroutine neighbor_lists() bind(C) end subroutine neighbor_lists end interface @@ -68,121 +78,4 @@ contains end subroutine find_cell -!=============================================================================== -! DISTANCE_TO_BOUNDARY calculates the distance to the nearest boundary for a -! particle 'p' traveling in a certain direction. For a cell in a subuniverse -! that has a parent cell, also include the surfaces of the edge of the universe. -!=============================================================================== - - subroutine distance_to_boundary(p, dist, surface_crossed, lattice_translation, & - next_level) - type(Particle), intent(inout) :: p - real(8), intent(out) :: dist - integer, intent(out) :: surface_crossed - integer, intent(out) :: lattice_translation(3) - integer, intent(out) :: next_level - - integer :: j - integer :: i_xyz(3) ! lattice indices - integer :: level_surf_cross ! surface crossed on current level - integer :: level_lat_trans(3) ! lattice translation on current level - real(8) :: xyz_t(3) ! local particle coordinates - real(8) :: d_lat ! distance to lattice boundary - real(8) :: d_surf ! distance to surface - real(8) :: xyz_cross(3) ! coordinates at projected surface crossing - real(8) :: surf_uvw(3) ! surface normal direction - type(Cell), pointer :: c - class(Lattice), pointer :: lat - - ! inialize distance to infinity (huge) - dist = INFINITY - d_lat = INFINITY - d_surf = INFINITY - lattice_translation(:) = [0, 0, 0] - - next_level = 0 - - ! Loop over each universe level - LEVEL_LOOP: do j = 1, p % n_coord - - ! get pointer to cell on this level - c => cells(p % coord(j) % cell) - - ! ======================================================================= - ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL - - call c % distance(p % coord(j) % xyz, p % coord(j) % uvw, p % surface, & - d_surf, level_surf_cross) - - ! ======================================================================= - ! FIND MINIMUM DISTANCE TO LATTICE SURFACES - - LAT_COORD: if (p % coord(j) % lattice /= NONE) then - lat => lattices(p % coord(j) % lattice) % obj - - i_xyz(1) = p % coord(j) % lattice_x - i_xyz(2) = p % coord(j) % lattice_y - i_xyz(3) = p % coord(j) % lattice_z - - LAT_TYPE: select type(lat) - - type is (RectLattice) - call lat % distance(p % coord(j) % xyz, p % coord(j) % uvw, & - i_xyz, d_lat, level_lat_trans) - - type is (HexLattice) LAT_TYPE - xyz_t(1) = p % coord(j-1) % xyz(1) - xyz_t(2) = p % coord(j-1) % xyz(2) - xyz_t(3) = p % coord(j) % xyz(3) - call lat % distance(xyz_t, p % coord(j) % uvw, & - i_xyz, d_lat, level_lat_trans) - end select LAT_TYPE - - if (d_lat < ZERO) then - call particle_mark_as_lost(p, "Particle " // trim(to_str(p % id)) & - //" had a negative distance to a lattice boundary. d = " & - //trim(to_str(d_lat))) - end if - end if LAT_COORD - - ! If the boundary on this lattice level is coincident with a boundary on - ! a higher level then we need to make sure that the higher level boundary - ! is selected. This logic must include consideration of floating point - ! precision. - if (d_surf < d_lat) then - if ((dist - d_surf)/dist >= FP_REL_PRECISION) then - dist = d_surf - - ! If the cell is not simple, it is possible that both the negative and - ! positive half-space were given in the region specification. Thus, we - ! have to explicitly check which half-space the particle would be - ! traveling into if the surface is crossed - if (.not. c % simple()) then - xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw - call surfaces(abs(level_surf_cross)) % normal(xyz_cross, surf_uvw) - if (dot_product(p % coord(j) % uvw, surf_uvw) > ZERO) then - surface_crossed = abs(level_surf_cross) - else - surface_crossed = -abs(level_surf_cross) - end if - else - surface_crossed = level_surf_cross - end if - - lattice_translation(:) = [0, 0, 0] - next_level = j - end if - else - if ((dist - d_lat)/dist >= FP_REL_PRECISION) then - dist = d_lat - surface_crossed = NONE - lattice_translation(:) = level_lat_trans - next_level = j - end if - end if - - end do LEVEL_LOOP - - end subroutine distance_to_boundary - end module geometry diff --git a/src/geometry.cpp b/src/geometry.cpp index 87bf0abfd..e7cbdd804 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -15,6 +15,8 @@ namespace openmc { std::vector overlap_check_count; +constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE + //============================================================================== extern "C" bool @@ -316,4 +318,101 @@ cross_lattice(Particle* p, int lattice_translation[3]) } } +//============================================================================== + +extern "C" void +distance_to_boundary(Particle* p, double* dist, int* surface_crossed, + int lattice_translation[3], int* next_level) +{ + *dist = INFINITY; + double d_lat = INFINITY; + double d_surf = INFINITY; + lattice_translation[0] = 0; + lattice_translation[1] = 0; + lattice_translation[2] = 0; + int32_t level_surf_cross; + std::array level_lat_trans; + + // Loop over each coordinate level. + for (int i = 0; i < p->n_coord; i++) { + Position r {p->coord[i].xyz}; + Direction u {p->coord[i].uvw}; + Cell& c {*global_cells[p->coord[i].cell-1]}; + + // Find the oncoming surface in this cell and the distance to it. + auto surface_distance = c.distance(r, u, p->surface); + d_surf = surface_distance.first; + level_surf_cross = surface_distance.second; + + // Find the distance to the next lattice tile crossing. + if (p->coord[i].lattice != F90_NONE) { + Lattice& lat {*lattices_c[p->coord[i].lattice-1]}; + std::array i_xyz {p->coord[i].lattice_x, p->coord[i].lattice_y, + p->coord[i].lattice_z}; + //TODO: refactor so both lattice use the same position argument (which + //also means the lat.type attribute can be removed) + std::pair> lattice_distance; + switch (lat.type) { + case LatticeType::rect: + lattice_distance = lat.distance(r, u, i_xyz); + break; + case LatticeType::hex: + Position r_hex {p->coord[i-1].xyz[0], p->coord[i-1].xyz[1], + p->coord[i].xyz[2]}; + lattice_distance = lat.distance(r_hex, u, i_xyz); + break; + } + d_lat = lattice_distance.first; + level_lat_trans = lattice_distance.second; + + if (d_lat < 0) { + std::stringstream err_msg; + err_msg << "Particle " << p->id + << " had a negative distance to a lattice boundary"; + p->mark_as_lost(err_msg); + } + } + + // If the boundary on this coordinate level is coincident with a boundary on + // a higher level then we need to make sure that the higher level boundary + // is selected. This logic must consider floating point precision. + if (d_surf < d_lat) { + if (*dist == INFINITY || ((*dist) - d_surf)/(*dist) >= FP_REL_PRECISION) { + *dist = d_surf; + + // If the cell is not simple, it is possible that both the negative and + // positive half-space were given in the region specification. Thus, we + // have to explicitly check which half-space the particle would be + // traveling into if the surface is crossed + if (c.simple) { + *surface_crossed = level_surf_cross; + } else { + Position r_hit = r + d_surf * u; + Surface& surf {*global_surfaces[std::abs(level_surf_cross)-1]}; + Direction norm = surf.normal(r_hit); + if (u.dot(norm) > 0) { + *surface_crossed = std::abs(level_surf_cross); + } else { + *surface_crossed = -std::abs(level_surf_cross); + } + } + + lattice_translation[0] = 0; + lattice_translation[1] = 0; + lattice_translation[2] = 0; + *next_level = i + 1; + } + } else { + if (*dist == INFINITY || ((*dist) - d_lat)/(*dist) >= FP_REL_PRECISION) { + *dist = d_lat; + *surface_crossed = F90_NONE; + lattice_translation[0] = level_lat_trans[0]; + lattice_translation[1] = level_lat_trans[1]; + lattice_translation[2] = level_lat_trans[2]; + *next_level = i + 1; + } + } + } +} + } // namespace openmc diff --git a/src/geometry.h b/src/geometry.h index 9a13ae220..742818579 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -33,6 +33,14 @@ find_cell(Particle* p, int search_surf); extern "C" void cross_lattice(Particle* p, int lattice_translation[3]); +//============================================================================== +//! Find the next boundary a particle will intersect. +//============================================================================== + +extern "C" void +distance_to_boundary(Particle* p, double* dist, int* surface_crossed, + int lattice_translation[3], int* next_level); + } // namespace openmc #endif // OPENMC_GEOMETRY_H diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 0223a63ff..504534d30 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -97,23 +97,6 @@ module geometry_header real(C_DOUBLE) :: sqrtkT end function cell_sqrtkT_c - function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple) - import C_PTR, C_BOOL - type(C_PTR), intent(in), value :: cell_ptr - logical(C_BOOL) :: simple - end function cell_simple_c - - subroutine cell_distance_c(cell_ptr, xyz, uvw, on_surface, min_dist, & - i_surf) bind(C, name="cell_distance") - import C_PTR, C_INT32_T, C_DOUBLE - type(C_PTR), intent(in), value :: cell_ptr - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - integer(C_INT32_T), intent(in), value :: on_surface - real(C_DOUBLE), intent(out) :: min_dist - integer(C_INT32_T), intent(out) :: i_surf - end subroutine cell_distance_c - function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") & result(offset) import C_PTR, C_INT, C_INT32_T @@ -148,17 +131,6 @@ module geometry_header logical(C_BOOL) :: is_valid end function lattice_are_valid_indices_c - subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) & - bind(C, name='lattice_distance') - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), intent(in), value :: lat_ptr - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - integer(C_INT), intent(in) :: i_xyz(3) - real(C_DOUBLE), intent(out) :: d - integer(C_INT), intent(out) :: lattice_trans(3) - end subroutine lattice_distance_c - subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) & bind(C, name='lattice_get_indices') import C_PTR, C_INT, C_DOUBLE @@ -223,7 +195,6 @@ module geometry_header contains procedure :: id => lattice_id procedure :: are_valid_indices => lattice_are_valid_indices - procedure :: distance => lattice_distance procedure :: get => lattice_get procedure :: get_indices => lattice_get_indices procedure :: get_local_xyz => lattice_get_local_xyz @@ -277,8 +248,6 @@ module geometry_header procedure :: material => cell_material procedure :: sqrtkT_size => cell_sqrtkT_size procedure :: sqrtkT => cell_sqrtkT - procedure :: simple => cell_simple - procedure :: distance => cell_distance procedure :: offset => cell_offset procedure :: to_hdf5 => cell_to_hdf5 @@ -314,16 +283,6 @@ contains is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) end function lattice_are_valid_indices - subroutine lattice_distance(this, xyz, uvw, i_xyz, d, lattice_trans) - class(Lattice), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - integer(C_INT), intent(in) :: i_xyz(3) - real(C_DOUBLE), intent(out) :: d - integer(C_INT), intent(out) :: lattice_trans(3) - call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans) - end subroutine lattice_distance - function lattice_get(this, i_xyz) result(univ) class(Lattice), intent(in) :: this integer(C_INT), intent(in) :: i_xyz(3) @@ -431,22 +390,6 @@ contains sqrtkT = cell_sqrtkT_c(this % ptr, i) end function cell_sqrtkT - function cell_simple(this) result(simple) - class(Cell), intent(in) :: this - logical(C_BOOL) :: simple - simple = cell_simple_c(this % ptr) - end function cell_simple - - subroutine cell_distance(this, xyz, uvw, on_surface, min_dist, i_surf) - class(Cell), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - integer(C_INT32_T), intent(in) :: on_surface - real(C_DOUBLE), intent(out) :: min_dist - integer(C_INT32_T), intent(out) :: i_surf - call cell_distance_c(this % ptr, xyz, uvw, on_surface, min_dist, i_surf) - end subroutine cell_distance - function cell_offset(this, map) result(offset) class(Cell), intent(in) :: this integer(C_INT), intent(in) :: map diff --git a/src/lattice.cpp b/src/lattice.cpp index 8797d0fba..46f88e351 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -137,6 +137,8 @@ Lattice::to_hdf5(hid_t lattices_group) const RectLattice::RectLattice(pugi::xml_node lat_node) : Lattice {lat_node} { + type = LatticeType::rect; + // Read the number of lattice cells in each dimension. std::string dimension_str {get_node_value(lat_node, "dimension")}; std::vector dimension_words {split(dimension_str)}; @@ -399,6 +401,8 @@ RectLattice::to_hdf5_inner(hid_t lat_group) const HexLattice::HexLattice(pugi::xml_node lat_node) : Lattice {lat_node} { + type = LatticeType::hex; + // Read the number of lattice cells in each dimension. n_rings = std::stoi(get_node_value(lat_node, "n_rings")); if (check_for_node(lat_node, "n_axial")) { @@ -889,18 +893,6 @@ extern "C" { bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) {return lat->are_valid_indices(i_xyz);} - void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], - const int i_xyz[3], double *d, int lattice_trans[3]) - { - Position r {xyz}; - Direction u {uvw}; - std::pair> ld {lat->distance(r, u, i_xyz)}; - *d = ld.first; - lattice_trans[0] = ld.second[0]; - lattice_trans[1] = ld.second[1]; - lattice_trans[2] = ld.second[2]; - } - void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) { Position r {xyz}; diff --git a/src/lattice.h b/src/lattice.h index df7c699d8..758acfa82 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -22,6 +22,10 @@ namespace openmc { constexpr int32_t NO_OUTER_UNIVERSE{-1}; +enum class LatticeType { + rect, hex +}; + //============================================================================== // Global variables //============================================================================== @@ -44,6 +48,7 @@ class Lattice public: int32_t id; //!< Universe ID number std::string name; //!< User-defined name + LatticeType type; std::vector universes; //!< Universes filling each lattice tile int32_t outer {NO_OUTER_UNIVERSE}; //!< Universe tiled outside the lattice std::vector offsets; //!< Distribcell offset table @@ -100,6 +105,13 @@ public: distance(Position r, Direction u, const int i_xyz[3]) const = 0; + std::pair> + distance(Position r, Direction u, std::array i_xyz) const + { + int i_xyz_[3] {i_xyz[0], i_xyz[1], i_xyz[2]}; + return distance(r, u, i_xyz_); + } + //! \brief Find the lattice tile indices for a given point. //! \param r A 3D Cartesian coordinate. //! \return An array containing the indices of a lattice tile.