crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move distance_to_boundary to C++

Browse source
This commit is contained in:
Sterling Harper 2018-08-21 11:02:14 -04:00
parent bc1c1e4916
commit 8bb17a481e
7 changed files with 133 additions and 198 deletions
Download patch file Download diff file Expand all files Collapse all files

12
src/cell.cpp
View file

@ -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<double, int32_t> 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);}

127
src/geometry.F90
View file

@ -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

99
src/geometry.cpp
View file

@ -15,6 +15,8 @@ namespace openmc {
std::vector<int64_t> 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<int, 3> 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<int, 3> 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<double, std::array<int, 3>> 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

8
src/geometry.h
View file

@ -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

57
src/geometry_header.F90
View file

@ -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

16
src/lattice.cpp
View file

@ -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<std::string> 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<double, std::array<int, 3>> 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};

12
src/lattice.h
View file

@ -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<int32_t> universes; //!< Universes filling each lattice tile
int32_t outer {NO_OUTER_UNIVERSE}; //!< Universe tiled outside the lattice
std::vector<int32_t> offsets; //!< Distribcell offset table
@ -100,6 +105,13 @@ public:
distance(Position r, Direction u, const int i_xyz[3]) const
= 0;
std::pair<double, std::array<int, 3>>
distance(Position r, Direction u, std::array<int, 3> 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.