Have distance_to_boundary return a struct

This commit is contained in:
Paul Romano 2019-03-21 07:43:42 -05:00
parent 2e94189b1b
commit c60ecd7e8f
3 changed files with 52 additions and 54 deletions

View file

@ -1,6 +1,8 @@
#ifndef OPENMC_GEOMETRY_H
#define OPENMC_GEOMETRY_H
#include <array>
#include <cmath>
#include <cstdint>
#include <vector>
@ -22,12 +24,22 @@ extern std::vector<int64_t> overlap_check_count;
} // namespace model
//==============================================================================
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface_index {0}; //!< if boundary is surface, index in surfaces vector
int coord_level; //!< coordinate level after crossing boundary
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p);
bool check_cell_overlap(Particle* p);
//==============================================================================
//! Locate a particle in the geometry tree and set its geometry data fields.
@ -41,23 +53,19 @@ check_cell_overlap(Particle* p);
//! valid geometry coordinate stack.
//==============================================================================
extern "C" bool
find_cell(Particle* p, bool use_neighbor_lists);
bool find_cell(Particle* p, bool use_neighbor_lists);
//==============================================================================
//! Move a particle into a new lattice tile.
//==============================================================================
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3]);
void cross_lattice(Particle* p, const BoundaryInfo& boundary);
//==============================================================================
//! 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);
BoundaryInfo distance_to_boundary(Particle* p);
} // namespace openmc

View file

@ -32,8 +32,7 @@ std::vector<int64_t> overlap_check_count;
// Non-member functions
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p)
bool check_cell_overlap(Particle* p)
{
int n_coord = p->n_coord_;
@ -246,7 +245,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//==============================================================================
extern "C" bool
bool
find_cell(Particle* p, bool use_neighbor_lists)
{
// Determine universe (if not yet set, use root universe).
@ -288,8 +287,8 @@ find_cell(Particle* p, bool use_neighbor_lists)
//==============================================================================
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3])
void
cross_lattice(Particle* p, const BoundaryInfo& boundary)
{
auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]};
@ -303,9 +302,9 @@ cross_lattice(Particle* p, int lattice_translation[3])
}
// Set the lattice indices.
p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2];
p->coord_[p->n_coord_-1].lattice_x += boundary.lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += boundary.lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += boundary.lattice_translation[2];
std::array<int, 3> i_xyz {p->coord_[p->n_coord_-1].lattice_x,
p->coord_[p->n_coord_-1].lattice_y,
p->coord_[p->n_coord_-1].lattice_z};
@ -346,16 +345,11 @@ 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)
BoundaryInfo distance_to_boundary(Particle* p)
{
*dist = INFINITY;
BoundaryInfo info;
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;
@ -402,43 +396,43 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
// 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.
double& d = info.distance;
if (d_surf < d_lat) {
if (*dist == INFINITY || ((*dist) - d_surf)/(*dist) >= FP_REL_PRECISION) {
*dist = d_surf;
if (d == INFINITY || (d - d_surf)/d >= FP_REL_PRECISION) {
d = 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;
info.surface_index = level_surf_cross;
} else {
Position r_hit = r + d_surf * u;
Surface& surf {*model::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);
info.surface_index = std::abs(level_surf_cross);
} else {
*surface_crossed = -std::abs(level_surf_cross);
info.surface_index = -std::abs(level_surf_cross);
}
}
lattice_translation[0] = 0;
lattice_translation[1] = 0;
lattice_translation[2] = 0;
*next_level = i + 1;
info.lattice_translation[0] = 0;
info.lattice_translation[1] = 0;
info.lattice_translation[2] = 0;
info.coord_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;
if (d == INFINITY || (d - d_lat)/d >= FP_REL_PRECISION) {
d = d_lat;
info.surface_index = 0;
info.lattice_translation = level_lat_trans;
info.coord_level = i + 1;
}
}
}
return info;
}
//==============================================================================

View file

@ -202,12 +202,7 @@ Particle::transport()
}
// Find the distance to the nearest boundary
double d_boundary;
int surface_crossed;
int lattice_translation[3];
int next_level;
distance_to_boundary(this, &d_boundary, &surface_crossed,
lattice_translation, &next_level);
auto boundary = distance_to_boundary(this);
// Sample a distance to collision
double d_collision;
@ -221,7 +216,7 @@ Particle::transport()
}
// Select smaller of the two distances
double distance = std::min(d_boundary, d_collision);
double distance = std::min(boundary.distance, d_collision);
// Advance particle
for (int j = 0; j < n_coord_; ++j) {
@ -244,11 +239,13 @@ Particle::transport()
score_track_derivative(this, distance);
}
if (d_collision > d_boundary) {
if (d_collision > boundary.distance) {
// ====================================================================
// PARTICLE CROSSES SURFACE
if (next_level > 0) n_coord_ = next_level;
// Set surface that particle is on and adjust coordinate levels
surface_ = boundary.surface_index;
n_coord_ = boundary.coord_level;
// Saving previous cell data
for (int j = 0; j < n_coord_; ++j) {
@ -256,15 +253,14 @@ Particle::transport()
}
n_coord_last_ = n_coord_;
if (lattice_translation[0] != 0 || lattice_translation[1] != 0 ||
lattice_translation[2] != 0) {
if (boundary.lattice_translation[0] != 0 ||
boundary.lattice_translation[1] != 0 ||
boundary.lattice_translation[2] != 0) {
// Particle crosses lattice boundary
surface_ = 0;
cross_lattice(this, lattice_translation);
cross_lattice(this, boundary);
event_ = EVENT_LATTICE;
} else {
// Particle crosses surface
surface_ = surface_crossed;
this->cross_surface();
event_ = EVENT_SURFACE;
}