Boundary info accessors (#3496)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Patrick Shriwise 2025-07-18 02:51:30 -05:00 committed by GitHub
parent 4483583ddd
commit 637e04a9ba
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
9 changed files with 81 additions and 63 deletions

View file

@ -210,23 +210,41 @@ struct CacheDataMG {
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface {
SURFACE_NONE}; //!< surface token, non-zero if boundary is surface
int coord_level; //!< coordinate level after crossing boundary
array<int, 3>
lattice_translation {}; //!< which way lattice indices will change
class BoundaryInfo {
public:
void reset()
{
distance = INFINITY;
surface = SURFACE_NONE;
coord_level = 0;
lattice_translation = {0, 0, 0};
distance_ = INFINITY;
surface_ = SURFACE_NONE;
coord_level_ = 0;
lattice_translation_ = {0, 0, 0};
}
double& distance() { return distance_; }
const double& distance() const { return distance_; }
int& surface() { return surface_; }
const int& surface() const { return surface_; }
int coord_level() const { return coord_level_; }
int& coord_level() { return coord_level_; }
array<int, 3>& lattice_translation() { return lattice_translation_; }
const array<int, 3>& lattice_translation() const
{
return lattice_translation_;
}
// TODO: off-by-one
int surface_index() const { return std::abs(surface) - 1; }
int surface_index() const { return std::abs(surface()) - 1; }
private:
// Data members
double distance_ {INFINITY}; //!< distance to nearest boundary
int surface_ {
SURFACE_NONE}; //!< surface token, non-zero if boundary is surface
int coord_level_ {0}; //!< coordinate level after crossing boundary
array<int, 3> lattice_translation_ {
0, 0, 0}; //!< which way lattice indices will change
};
/*

View file

@ -114,7 +114,7 @@ void process_advance_particle_events()
p.event_advance();
if (!p.alive())
continue;
if (p.collision_distance() > p.boundary().distance) {
if (p.collision_distance() > p.boundary().distance()) {
simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
} else {
simulation::collision_queue.thread_safe_append({p, buffer_idx});

View file

@ -310,9 +310,9 @@ void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
}
// Set the lattice indices.
coord.lattice_index()[0] += boundary.lattice_translation[0];
coord.lattice_index()[1] += boundary.lattice_translation[1];
coord.lattice_index()[2] += boundary.lattice_translation[2];
coord.lattice_index()[0] += boundary.lattice_translation()[0];
coord.lattice_index()[1] += boundary.lattice_translation()[1];
coord.lattice_index()[2] += boundary.lattice_translation()[2];
// Set the new coordinate position.
const auto& upper_coord {p.coord(p.n_coord() - 2)};
@ -410,7 +410,7 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
// 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;
double& d = info.distance();
if (d_surf < d_lat - FP_COINCIDENT) {
if (d == INFINITY || (d - d_surf) / d >= FP_REL_PRECISION) {
// Update closest distance
@ -421,29 +421,29 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
// have to explicitly check which half-space the particle would be
// traveling into if the surface is crossed
if (c.is_simple() || d == INFTY) {
info.surface = level_surf_cross;
info.surface() = 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) {
info.surface = std::abs(level_surf_cross);
info.surface() = std::abs(level_surf_cross);
} else {
info.surface = -std::abs(level_surf_cross);
info.surface() = -std::abs(level_surf_cross);
}
}
info.lattice_translation[0] = 0;
info.lattice_translation[1] = 0;
info.lattice_translation[2] = 0;
info.coord_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 (d == INFINITY || (d - d_lat) / d >= FP_REL_PRECISION) {
d = d_lat;
info.surface = SURFACE_NONE;
info.lattice_translation = level_lat_trans;
info.coord_level = i + 1;
info.surface() = SURFACE_NONE;
info.lattice_translation() = level_lat_trans;
info.coord_level() = i + 1;
}
}
}

View file

@ -385,7 +385,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
BoundaryInfo boundary = distance_to_boundary(p);
// Advance particle forward
double distance = std::min(boundary.distance, max_distance);
double distance = std::min(boundary.distance(), max_distance);
p.move_distance(distance);
// Determine what mesh elements were crossed by particle
@ -416,12 +416,12 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
p.n_coord_last() = p.n_coord();
// Set surface that particle is on and adjust coordinate levels
p.surface() = boundary.surface;
p.n_coord() = boundary.coord_level;
p.surface() = boundary.surface();
p.n_coord() = boundary.coord_level();
if (boundary.lattice_translation[0] != 0 ||
boundary.lattice_translation[1] != 0 ||
boundary.lattice_translation[2] != 0) {
if (boundary.lattice_translation()[0] != 0 ||
boundary.lattice_translation()[1] != 0 ||
boundary.lattice_translation()[2] != 0) {
// Particle crosses lattice boundary
cross_lattice(p, boundary);
} else {

View file

@ -242,7 +242,7 @@ void Particle::event_advance()
}
// Select smaller of the two distances
double distance = std::min(boundary().distance, collision_distance());
double distance = std::min(boundary().distance(), collision_distance());
// Advance particle in space and time
// Short-term solution until the surface source is revised and we can use
@ -298,12 +298,12 @@ void Particle::event_cross_surface()
n_coord_last() = n_coord();
// Set surface that particle is on and adjust coordinate levels
surface() = boundary().surface;
n_coord() = boundary().coord_level;
surface() = boundary().surface();
n_coord() = boundary().coord_level();
if (boundary().lattice_translation[0] != 0 ||
boundary().lattice_translation[1] != 0 ||
boundary().lattice_translation[2] != 0) {
if (boundary().lattice_translation()[0] != 0 ||
boundary().lattice_translation()[1] != 0 ||
boundary().lattice_translation()[2] != 0) {
// Particle crosses lattice boundary
bool verbose = settings::verbosity >= 10 || trace();

View file

@ -65,22 +65,22 @@ void GeometryState::advance_to_boundary_from_void()
for (auto c_i : root_universe->cells_) {
auto dist =
model::cells.at(c_i)->distance(root_coord.r(), root_coord.u(), 0, this);
if (dist.first < boundary().distance) {
boundary().distance = dist.first;
boundary().surface = dist.second;
if (dist.first < boundary().distance()) {
boundary().distance() = dist.first;
boundary().surface() = dist.second;
}
}
// if no intersection or near-infinite intersection, reset
// boundary information
if (boundary().distance > 1e300) {
boundary().distance = INFTY;
boundary().surface = SURFACE_NONE;
if (boundary().distance() > 1e300) {
boundary().distance() = INFTY;
boundary().surface() = SURFACE_NONE;
return;
}
// move the particle up to (and just past) the boundary
move_distance(boundary().distance + TINY_BIT);
move_distance(boundary().distance() + TINY_BIT);
}
void GeometryState::move_distance(double length)

View file

@ -1634,7 +1634,7 @@ void Ray::trace()
break;
// if there is no intersection with the model, we're done
if (boundary().surface == SURFACE_NONE)
if (boundary().surface() == SURFACE_NONE)
return;
event_counter_++;
@ -1646,10 +1646,10 @@ void Ray::trace()
// Call the specialized logic for this type of ray. This is for the
// intersection for the first intersection if we had one.
if (boundary().surface != SURFACE_NONE) {
if (boundary().surface() != SURFACE_NONE) {
// set the geometry state's surface attribute to be used for
// surface normal computation
surface() = boundary().surface;
surface() = boundary().surface();
on_intersection();
if (stop_)
return;
@ -1674,37 +1674,37 @@ void Ray::trace()
// if we hit the edge of the model, so stop
// the particle in that case. Also, just exit
// if a negative distance was somehow computed.
if (boundary().distance == INFTY || boundary().distance == INFINITY ||
boundary().distance < 0) {
if (boundary().distance() == INFTY || boundary().distance() == INFINITY ||
boundary().distance() < 0) {
return;
}
// See below comment where call_on_intersection is checked in an
// if statement for an explanation of this.
bool call_on_intersection {true};
if (boundary().distance < 10 * TINY_BIT) {
if (boundary().distance() < 10 * TINY_BIT) {
call_on_intersection = false;
}
// DAGMC surfaces expect us to go a little bit further than the advance
// distance to properly check cell inclusion.
boundary().distance += TINY_BIT;
boundary().distance() += TINY_BIT;
// Advance particle, prepare for next intersection
for (int lev = 0; lev < n_coord(); ++lev) {
coord(lev).r() += boundary().distance * coord(lev).u();
coord(lev).r() += boundary().distance() * coord(lev).u();
}
surface() = boundary().surface;
surface() = boundary().surface();
n_coord_last() = n_coord();
n_coord() = boundary().coord_level;
if (boundary().lattice_translation[0] != 0 ||
boundary().lattice_translation[1] != 0 ||
boundary().lattice_translation[2] != 0) {
n_coord() = boundary().coord_level();
if (boundary().lattice_translation()[0] != 0 ||
boundary().lattice_translation()[1] != 0 ||
boundary().lattice_translation()[2] != 0) {
cross_lattice(*this, boundary(), settings::verbosity >= 10);
}
// Record how far the ray has traveled
traversal_distance_ += boundary().distance;
traversal_distance_ += boundary().distance();
inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
// Call the specialized logic for this type of ray. Note that we do not

View file

@ -281,7 +281,7 @@ void RandomRay::event_advance_ray()
{
// Find the distance to the nearest boundary
boundary() = distance_to_boundary(*this);
double distance = boundary().distance;
double distance = boundary().distance();
if (distance < 0.0) {
mark_as_lost("Negative transport distance detected for particle " +

View file

@ -790,7 +790,7 @@ void transport_history_based_single_particle(Particle& p)
p.event_advance();
}
if (p.alive()) {
if (p.collision_distance() > p.boundary().distance) {
if (p.collision_distance() > p.boundary().distance()) {
p.event_cross_surface();
} else if (p.alive()) {
p.event_collide();