mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 13:15:39 -04:00
Add C++ HexLattice distance_to_boundary
This commit is contained in:
Sterling Harper
parent
81cc98010b
commit
ec14970caf
4 changed files with 145 additions and 134 deletions
134
src/geometry.F90
134
src/geometry.F90
|
|
@ -372,17 +372,9 @@ contains
|
|||
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) :: x,y,z ! particle coordinates
|
||||
real(8) :: xyz_t(3) ! local particle coordinates
|
||||
real(8) :: beta, gama ! skewed particle coordiantes
|
||||
real(8) :: u,v,w ! particle directions
|
||||
real(8) :: beta_dir ! skewed particle direction
|
||||
real(8) :: gama_dir ! skewed particle direction
|
||||
real(8) :: edge ! distance to oncoming edge
|
||||
real(8) :: d ! evaluated distance
|
||||
real(8) :: d_lat ! distance to lattice boundary
|
||||
real(8) :: d_surf ! distance to surface
|
||||
real(8) :: x0,y0,z0 ! coefficients for surface
|
||||
real(8) :: xyz_cross(3) ! coordinates at projected surface crossing
|
||||
real(8) :: surf_uvw(3) ! surface normal direction
|
||||
type(Cell), pointer :: c
|
||||
|
|
@ -402,11 +394,6 @@ contains
|
|||
! get pointer to cell on this level
|
||||
c => cells(p % coord(j) % cell)
|
||||
|
||||
! copy directional cosines
|
||||
u = p % coord(j) % uvw(1)
|
||||
v = p % coord(j) % uvw(2)
|
||||
w = p % coord(j) % uvw(3)
|
||||
|
||||
! =======================================================================
|
||||
! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL
|
||||
|
||||
|
|
@ -419,123 +406,22 @@ contains
|
|||
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, &
|
||||
d_lat, level_lat_trans)
|
||||
i_xyz, d_lat, level_lat_trans)
|
||||
|
||||
type is (HexLattice) LAT_TYPE
|
||||
! Copy local coordinates.
|
||||
z = p % coord(j) % xyz(3)
|
||||
i_xyz(1) = p % coord(j) % lattice_x
|
||||
i_xyz(2) = p % coord(j) % lattice_y
|
||||
i_xyz(3) = p % coord(j) % lattice_z
|
||||
|
||||
! Compute velocities along the hexagonal axes.
|
||||
beta_dir = u*sqrt(THREE)/TWO + v/TWO
|
||||
gama_dir = u*sqrt(THREE)/TWO - v/TWO
|
||||
|
||||
! Note that hexagonal lattice distance calculations are performed
|
||||
! using the particle's coordinates relative to the neighbor lattice
|
||||
! cells, not relative to the particle's current cell. This is done
|
||||
! because there is significant disagreement between neighboring cells
|
||||
! on where the lattice boundary is due to the worse finite precision
|
||||
! of hex lattices.
|
||||
|
||||
! Upper right and lower left sides.
|
||||
edge = -sign(lat % pitch(1)/TWO, beta_dir) ! Oncoming edge
|
||||
if (beta_dir > ZERO) then
|
||||
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, 0, 0])
|
||||
else
|
||||
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 0, 0])
|
||||
end if
|
||||
beta = xyz_t(1)*sqrt(THREE)/TWO + xyz_t(2)/TWO
|
||||
if (abs(beta - edge) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
else if (beta_dir == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (edge - beta)/beta_dir
|
||||
end if
|
||||
|
||||
d_lat = d
|
||||
if (beta_dir > 0) then
|
||||
level_lat_trans(:) = [1, 0, 0]
|
||||
else
|
||||
level_lat_trans(:) = [-1, 0, 0]
|
||||
end if
|
||||
|
||||
! Lower right and upper left sides.
|
||||
edge = -sign(lat % pitch(1)/TWO, gama_dir) ! Oncoming edge
|
||||
if (gama_dir > ZERO) then
|
||||
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, -1, 0])
|
||||
else
|
||||
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 1, 0])
|
||||
end if
|
||||
gama = xyz_t(1)*sqrt(THREE)/TWO - xyz_t(2)/TWO
|
||||
if (abs(gama - edge) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
else if (gama_dir == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (edge - gama)/gama_dir
|
||||
end if
|
||||
|
||||
if (d < d_lat) then
|
||||
d_lat = d
|
||||
if (gama_dir > 0) then
|
||||
level_lat_trans(:) = [1, -1, 0]
|
||||
else
|
||||
level_lat_trans(:) = [-1, 1, 0]
|
||||
end if
|
||||
end if
|
||||
|
||||
! Upper and lower sides.
|
||||
edge = -sign(lat % pitch(1)/TWO, v) ! Oncoming edge
|
||||
if (v > ZERO) then
|
||||
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, 1, 0])
|
||||
else
|
||||
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, -1, 0])
|
||||
end if
|
||||
if (abs(xyz_t(2) - edge) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
else if (v == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (edge - xyz_t(2))/v
|
||||
end if
|
||||
|
||||
if (d < d_lat) then
|
||||
d_lat = d
|
||||
if (v > 0) then
|
||||
level_lat_trans(:) = [0, 1, 0]
|
||||
else
|
||||
level_lat_trans(:) = [0, -1, 0]
|
||||
end if
|
||||
end if
|
||||
|
||||
! Top and bottom sides.
|
||||
if (lat % is_3d) then
|
||||
z0 = sign(lat % pitch(2) * HALF, w)
|
||||
|
||||
if (abs(z - z0) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
elseif (w == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (z0 - z)/w
|
||||
end if
|
||||
|
||||
if (d < d_lat) then
|
||||
d_lat = d
|
||||
if (w > 0) then
|
||||
level_lat_trans(:) = [0, 0, 1]
|
||||
else
|
||||
level_lat_trans(:) = [0, 0, -1]
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
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
|
||||
|
|
|
|||
|
|
@ -95,13 +95,14 @@ module geometry_header
|
|||
integer(C_INT32_T) :: id
|
||||
end function lattice_id_c
|
||||
|
||||
subroutine lattice_distance_c(lat_ptr, xyz, uvw, d, lattice_trans) &
|
||||
subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) &
|
||||
bind(C, name='lattice_distance')
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
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
|
||||
|
|
@ -293,13 +294,14 @@ contains
|
|||
id = lattice_id_c(this % ptr)
|
||||
end function lattice_id
|
||||
|
||||
subroutine lattice_distance(this, xyz, uvw, d, lattice_trans)
|
||||
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, d, lattice_trans)
|
||||
call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans)
|
||||
end subroutine lattice_distance
|
||||
|
||||
subroutine lattice_to_hdf5(this, group)
|
||||
|
|
|
|||
127
src/lattice.cpp
127
src/lattice.cpp
|
|
@ -132,7 +132,8 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
|
|||
}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
RectLattice::distance(const double xyz[3], const double uvw[3]) const
|
||||
RectLattice::distance(const double xyz[3], const double uvw[3],
|
||||
const int i_xyz[3]) const
|
||||
{
|
||||
// Get short aliases to the coordinates.
|
||||
double x {xyz[0]};
|
||||
|
|
@ -342,8 +343,126 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
|
|||
}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
HexLattice::distance(const double xyz[3], const double uvw[3]) const
|
||||
HexLattice::distance(const double xyz[3], const double uvw[3],
|
||||
const int i_xyz[3]) const
|
||||
{
|
||||
// Compute the direction on the hexagonal basis.
|
||||
double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0;
|
||||
double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0;
|
||||
|
||||
// Note that hexagonal lattice distance calculations are performed
|
||||
// using the particle's coordinates relative to the neighbor lattice
|
||||
// cells, not relative to the particle's current cell. This is done
|
||||
// because there is significant disagreement between neighboring cells
|
||||
// on where the lattice boundary is due to finite precision issues.
|
||||
|
||||
// Upper-right and lower-left sides.
|
||||
double d {INFTY};
|
||||
std::array<int, 3> lattice_trans;
|
||||
double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge
|
||||
std::array<double, 3> xyz_t;
|
||||
if (beta_dir > 0) {
|
||||
const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
} else {
|
||||
const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
}
|
||||
double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0;
|
||||
if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
|
||||
d = (edge - beta) / beta_dir;
|
||||
if (beta_dir > 0) {
|
||||
lattice_trans = {1, 0, 0};
|
||||
} else {
|
||||
lattice_trans = {-1, 0, 0};
|
||||
}
|
||||
}
|
||||
|
||||
// Lower-right and upper-left sides.
|
||||
edge = -copysign(0.5*pitch[0], gamma_dir);
|
||||
if (gamma_dir > 0) {
|
||||
const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
} else {
|
||||
const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
}
|
||||
double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0;
|
||||
if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
|
||||
double this_d = (edge - gamma) / gamma_dir;
|
||||
if (this_d < d) {
|
||||
if (gamma_dir > 0) {
|
||||
lattice_trans = {1, -1, 0};
|
||||
} else {
|
||||
lattice_trans = {-1, 1, 0};
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
|
||||
// Upper and lower sides.
|
||||
edge = -copysign(0.5*pitch[0], uvw[1]);
|
||||
if (uvw[1] > 0) {
|
||||
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
} else {
|
||||
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
}
|
||||
if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) {
|
||||
double this_d = (edge - xyz_t[1]) / uvw[1];
|
||||
if (this_d < d) {
|
||||
if (uvw[1] > 0) {
|
||||
lattice_trans = {0, 1, 0};
|
||||
} else {
|
||||
lattice_trans = {0, -1, 0};
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
|
||||
// Top and bottom sides
|
||||
if (is_3d) {
|
||||
double z {xyz[2]};
|
||||
double w {uvw[2]};
|
||||
double z0 {copysign(0.5 * pitch[1], w)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
|
||||
double this_d = (z0 - z) / w;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (w > 0) {
|
||||
lattice_trans = {0, 0, 1};
|
||||
} else {
|
||||
lattice_trans = {0, 0, -1};
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return {d, lattice_trans};
|
||||
}
|
||||
|
||||
std::array<double, 3>
|
||||
HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
|
||||
{
|
||||
std::array<double, 3> local_xyz;
|
||||
|
||||
// x_l = x_g - (center + pitch_x*cos(30)*index_x)
|
||||
local_xyz[0] = global_xyz[0] - (center[0]
|
||||
+ std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings) * pitch[0]);
|
||||
// y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
|
||||
local_xyz[1] = global_xyz[1] - (center[1]
|
||||
+ (i_xyz[1] - n_rings) * pitch[0]
|
||||
+ (i_xyz[0] - n_rings) * pitch[0] / 2.0);
|
||||
if (is_3d) {
|
||||
local_xyz[2] = global_xyz[2] - center[2]
|
||||
+ (0.5 * n_axial - i_xyz[2] + 0.5) * pitch[1];
|
||||
} else {
|
||||
local_xyz[2] = global_xyz[2];
|
||||
}
|
||||
|
||||
return local_xyz;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -382,9 +501,9 @@ extern "C" {
|
|||
int32_t lattice_id(Lattice *lat) {return lat->id;}
|
||||
|
||||
void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3],
|
||||
double *d, int lattice_trans[3])
|
||||
const int i_xyz[3], double *d, int lattice_trans[3])
|
||||
{
|
||||
std::pair<double, std::array<int, 3>> ld {lat->distance(xyz, uvw)};
|
||||
std::pair<double, std::array<int, 3>> ld {lat->distance(xyz, uvw, i_xyz)};
|
||||
*d = ld.first;
|
||||
lattice_trans[0] = ld.second[0];
|
||||
lattice_trans[1] = ld.second[1];
|
||||
|
|
|
|||
|
|
@ -61,7 +61,8 @@ public:
|
|||
//virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
|
||||
|
||||
virtual std::pair<double, std::array<int, 3>>
|
||||
distance(const double xyz[3], const double uvw[3]) const = 0;
|
||||
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const
|
||||
= 0;
|
||||
|
||||
//virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0;
|
||||
|
||||
|
|
@ -87,7 +88,7 @@ public:
|
|||
virtual ~RectLattice() {}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
distance(const double xyz[3], const double uvw[3]) const;
|
||||
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
|
||||
|
||||
protected:
|
||||
std::array<int, 3> n_cells; //! Number of cells along each axis
|
||||
|
|
@ -103,7 +104,10 @@ public:
|
|||
virtual ~HexLattice() {}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
distance(const double xyz[3], const double uvw[3]) const;
|
||||
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
|
||||
|
||||
std::array<double, 3>
|
||||
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
|
||||
|
||||
protected:
|
||||
int n_rings; //! Number of radial tile positions
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue