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

Merge pull request #1200 from smharper/remove_mesh_search

Browse source 
Remove a search from tally mesh code
This commit is contained in:
Paul Romano 2019-03-19 12:48:18 -05:00 • committed by GitHub
commit 00343a327c
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 204 additions and 164 deletions
Download patch file Download diff file Expand all files Collapse all files

35
include/openmc/mesh.h
View file

@ -43,7 +43,7 @@ public:
// Methods
//! Determine which bins were crossed by a particle
//!
//
//! \param[in] p Particle to check
//! \param[out] bins Bins that were crossed
//! \param[out] lengths Fraction of tracklength in each bin
@ -51,50 +51,51 @@ public:
std::vector<double>& lengths) const;
//! Determine which surface bins were crossed by a particle
//!
//
//! \param[in] p Particle to check
//! \param[out] bins Surface bins that were crossed
void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
//! Get bin at a given position in space
//!
//
//! \param[in] r Position to get bin for
//! \return Mesh bin
int get_bin(Position r) const;
//! Get bin given mesh indices
//!
//
//! \param[in] Array of mesh indices
//! \return Mesh bin
int get_bin_from_indices(const int* ijk) const;
//! Get mesh indices given a position
//!
//
//! \param[in] r Position to get indices for
//! \param[out] ijk Array of mesh indices
//! \param[out] in_mesh Whether position is in mesh
void get_indices(Position r, int* ijk, bool* in_mesh) const;
//! Get mesh indices corresponding to a mesh bin
//!
//
//! \param[in] bin Mesh bin
//! \param[out] ijk Mesh indices
void get_indices_from_bin(int bin, int* ijk) const;
//! Check if a line connected by two points intersects the mesh
//!
//! \param[in] r0 Starting position
//! Check where a line segment intersects the mesh and if it intersects at all
//
//! \param[in,out] r0 In: starting position, out: intersection point
//! \param[in] r1 Ending position
//! \return Whether line connecting r0 and r1 intersects mesh
bool intersects(Position r0, Position r1) const;
//! \param[out] ijk Indices of the mesh bin containing the intersection point
//! \return Whether the line segment connecting r0 and r1 intersects mesh
bool intersects(Position& r0, Position r1, int* ijk) const;
//! Write mesh data to an HDF5 group
//!
//
//! \param[in] group HDF5 group
void to_hdf5(hid_t group) const;
//! Count number of bank sites in each mesh bin / energy bin
//!
//
//! \param[in] n Number of bank sites
//! \param[in] bank Array of bank sites
//! \param[in] n_energy Number of energies
@ -113,9 +114,9 @@ public:
xt::xarray<double> width_; //!< Width of each mesh element
private:
bool intersects_1d(Position r0, Position r1) const;
bool intersects_2d(Position r0, Position r1) const;
bool intersects_3d(Position r0, Position r1) const;
bool intersects_1d(Position& r0, Position r1, int* ijk) const;
bool intersects_2d(Position& r0, Position r1, int* ijk) const;
bool intersects_3d(Position& r0, Position r1, int* ijk) const;
};
//==============================================================================
@ -123,10 +124,12 @@ private:
//==============================================================================
//! Read meshes from either settings/tallies
//
//! \param[in] root XML node
void read_meshes(pugi::xml_node root);
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
void meshes_to_hdf5(hid_t group);

333
src/mesh.cpp
View file

@ -36,6 +36,32 @@ std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
//==============================================================================
// Helper functions
//==============================================================================
//! Update an intersection point if the given candidate is closer.
//
//! The first 6 arguments are coordinates for the starting point of a particle
//! and its intersection with a mesh surface. If the distance between these
//! two points is shorter than the given `min_distance`, then the `r` argument
//! will be updated to match the intersection point, and `min_distance` will
//! also be updated.
inline bool check_intersection_point(double x1, double x0, double y1,
double y0, double z1, double z0, Position& r, double& min_distance)
{
double dist = std::pow(x1-x0, 2) + std::pow(y1-y0, 2) + std::pow(z1-z0, 2);
if (dist < min_distance) {
r.x = x1;
r.y = y1;
r.z = z1;
min_distance = dist;
return true;
}
return false;
}
//==============================================================================
// RegularMesh implementation
//==============================================================================
@ -167,14 +193,14 @@ int RegularMesh::get_bin(Position r) const
int RegularMesh::get_bin_from_indices(const int* ijk) const
{
switch (n_dimension_) {
case 1:
return ijk[0] - 1;
case 2:
return (ijk[1] - 1)*shape_[0] + ijk[0] - 1;
case 3:
return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
default:
throw std::runtime_error{"Invalid number of mesh dimensions"};
case 1:
return ijk[0] - 1;
case 2:
return (ijk[1] - 1)*shape_[0] + ijk[0] - 1;
case 3:
return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
default:
throw std::runtime_error{"Invalid number of mesh dimensions"};
}
}
@ -204,45 +230,72 @@ void RegularMesh::get_indices_from_bin(int bin, int* ijk) const
}
}
bool RegularMesh::intersects(Position r0, Position r1) const
bool RegularMesh::intersects(Position& r0, Position r1, int* ijk) const
{
switch(n_dimension_) {
case 1:
return intersects_1d(r0, r1);
case 2:
return intersects_2d(r0, r1);
case 3:
return intersects_3d(r0, r1);
default:
throw std::runtime_error{"Invalid number of mesh dimensions."};
case 1:
return intersects_1d(r0, r1, ijk);
case 2:
return intersects_2d(r0, r1, ijk);
case 3:
return intersects_3d(r0, r1, ijk);
default:
throw std::runtime_error{"Invalid number of mesh dimensions."};
}
}
bool RegularMesh::intersects_1d(Position r0, Position r1) const
{
// Copy coordinates of mesh lower_left and upper_right
double left = lower_left_[0];
double right = upper_right_[0];
// Check if line intersects either left or right surface
if (r0.x < left) {
return r1.x > left;
} else if (r0.x < right) {
return r1.x < left || r1.x > right;
} else {
return r1.x < right;
}
}
bool RegularMesh::intersects_2d(Position r0, Position r1) const
bool RegularMesh::intersects_1d(Position& r0, Position r1, int* ijk) const
{
// Copy coordinates of starting point
double x0 = r0.x;
double y0 = r0.y;
double z0 = r0.z;
// Copy coordinates of ending point
double x1 = r1.x;
double y1 = r1.y;
double z1 = r1.z;
// Copy coordinates of mesh lower_left and upper_right
double xm0 = lower_left_[0];
double xm1 = upper_right_[0];
double min_dist = INFTY;
// Check if line intersects left surface -- calculate the intersection point
// (y,z)
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) {
ijk[0] = 1;
}
}
// Check if line intersects right surface -- calculate the intersection point
// (y,z)
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) {
ijk[0] = shape_[0];
}
}
return min_dist < INFTY;
}
bool RegularMesh::intersects_2d(Position& r0, Position r1, int* ijk) const
{
// Copy coordinates of starting point
double x0 = r0.x;
double y0 = r0.y;
double z0 = r0.z;
// Copy coordinates of ending point
double x1 = r1.x;
double y1 = r1.y;
double z1 = r1.z;
// Copy coordinates of mesh lower_left
double xm0 = lower_left_[0];
@ -252,44 +305,64 @@ bool RegularMesh::intersects_2d(Position r0, Position r1) const
double xm1 = upper_right_[0];
double ym1 = upper_right_[1];
// Check if line intersects left surface -- calculate the intersection point y
double min_dist = INFTY;
// Check if line intersects left surface -- calculate the intersection point
// (y,z)
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
if (yi >= ym0 && yi < ym1) {
return true;
if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) {
ijk[0] = 1;
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
}
}
}
// Check if line intersects back surface -- calculate the intersection point
// x
// (x,z)
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
if (xi >= xm0 && xi < xm1) {
return true;
if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) {
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
ijk[1] = 1;
}
}
}
// Check if line intersects right surface -- calculate the intersection
// point y
// Check if line intersects right surface -- calculate the intersection point
// (y,z)
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
if (yi >= ym0 && yi < ym1) {
return true;
if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) {
ijk[0] = shape_[0];
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
}
}
}
// Check if line intersects front surface -- calculate the intersection point
// x
// (x,z)
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
if (xi >= xm0 && xi < xm1) {
return true;
if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) {
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
ijk[1] = shape_[1];
}
}
}
return false;
return min_dist < INFTY;
}
bool RegularMesh::intersects_3d(Position r0, Position r1) const
bool RegularMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
{
// Copy coordinates of starting point
double x0 = r0.x;
@ -311,13 +384,19 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) const
double ym1 = upper_right_[1];
double zm1 = upper_right_[2];
double min_dist = INFTY;
// Check if line intersects left surface -- calculate the intersection point
// (y,z)
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
return true;
if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
ijk[0] = 1;
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
}
}
}
@ -327,7 +406,11 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) const
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
return true;
if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
ijk[1] = 1;
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
}
}
}
@ -337,7 +420,11 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) const
double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
return true;
if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
ijk[2] = 1;
}
}
}
@ -347,7 +434,11 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) const
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
return true;
if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
ijk[0] = shape_[0];
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
}
}
}
@ -357,7 +448,11 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) const
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
return true;
if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
ijk[1] = shape_[1];
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
}
}
}
@ -367,32 +462,38 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) const
double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
return true;
if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
ijk[2] = shape_[2];
}
}
}
return false;
return min_dist < INFTY;
}
void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& lengths) const
{
constexpr int MAX_SEARCH_ITER = 100;
// ========================================================================
// Determine if the track intersects the tally mesh.
// Determine where the track intersects the mesh and if it intersects at all.
// Copy the starting and ending coordinates of the particle. Offset these
// just a bit for the purposes of determining if there was an intersection
// in case the mesh surfaces coincide with lattice/geometric surfaces which
// might produce finite-precision errors.
// Copy the starting and ending coordinates of the particle.
Position last_r {p->r_last_};
Position r {p->r()};
Direction u {p->u()};
// Compute the length of the entire track.
double total_distance = (r - last_r).norm();
// While determining if this track intersects the mesh, offset the starting
// and ending coords by a bit. This avoid finite-precision errors that can
// occur when the mesh surfaces coincide with lattice or geometric surfaces.
Position r0 = last_r + TINY_BIT*u;
Position r1 = r - TINY_BIT*u;
// Determine indices for starting and ending location.
// Determine the mesh indices for the starting and ending coords.
int n = n_dimension_;
int ijk0[n], ijk1[n];
bool start_in_mesh;
@ -400,99 +501,32 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
bool end_in_mesh;
get_indices(r1, ijk1, &end_in_mesh);
// Check if the track intersects any part of the mesh.
if (!start_in_mesh && !end_in_mesh) {
if (!intersects(r0, r1)) return;
// Reset coordinates and check for a mesh intersection if necessary.
if (start_in_mesh) {
// The initial coords lie in the mesh, use those coords for tallying.
r0 = last_r;
} else {
// The initial coords do not lie in the mesh. Check to see if the particle
// eventually intersects the mesh and compute the relevant coords and
// indices.
if (!intersects(r0, r1, ijk0)) return;
}
// ========================================================================
// Figure out which mesh cell to tally.
// Copy the un-modified coordinates the particle direction.
r0 = last_r;
r1 = r;
// Compute the length of the entire track.
double total_distance = (r1 - r0).norm();
// We are looking for the first valid mesh bin. Check to see if the
// particle starts inside the mesh.
if (!start_in_mesh) {
double d[n];
// The particle does not start in the mesh. Note that we nudged the
// start and end coordinates by a TINY_BIT each so we will have
// difficulty resolving tracks that are less than 2*TINY_BIT in length.
// If the track is that short, it is also insignificant so we can
// safely ignore it in the tallies.
if (total_distance < 2*TINY_BIT) return;
// The particle does not start in the mesh so keep iterating the ijk0
// indices to cross the nearest mesh surface until we've found a valid
// bin. MAX_SEARCH_ITER prevents an infinite loop.
int search_iter = 0;
int j;
bool in_mesh = true;
for (int i = 0; i < n; ++i) {
if (ijk0[i] < 1 || ijk0[i] > shape_[i]) {
in_mesh = false;
break;
}
}
while (!in_mesh) {
if (search_iter == MAX_SEARCH_ITER) {
warning("Failed to find a mesh intersection on a tally mesh filter.");
return;
}
for (j = 0; j < n; ++j) {
if (std::fabs(u[j]) < FP_PRECISION) {
d[j] = INFTY;
} else if (u[j] > 0.0) {
double xyz_cross = lower_left_[j] + ijk0[j] * width_[j];
d[j] = (xyz_cross - r0[j]) / u[j];
} else {
double xyz_cross = lower_left_[j] + (ijk0[j] - 1) * width_[j];
d[j] = (xyz_cross - r0[j]) / u[j];
}
}
j = std::min_element(d, d+n) - d;
if (u[j] > 0.0) {
++ijk0[j];
} else {
--ijk0[j];
}
++search_iter;
in_mesh = true;
for (int i = 0; i < n; ++i) {
if (ijk0[i] < 1 || ijk0[i] > shape_[i]) {
in_mesh = false;
break;
}
}
}
// Advance position
r0 += d[j] * u;
}
// ========================================================================
// Find which mesh cells are traversed and the length of each traversal.
while (true) {
// ========================================================================
// Compute the length of the track segment in the each mesh cell and return
if (std::equal(ijk0, ijk0+n, ijk1)) {
// The track ends in this cell. Use the particle end location rather
// than the mesh surface.
// than the mesh surface and stop iterating.
double distance = (r1 - r0).norm();
bins.push_back(get_bin_from_indices(ijk0));
lengths.push_back(distance / total_distance);
break;
}
// The track exits this cell. Determine the distance to the closest mesh
// surface.
// The track exits this cell. Determine the distance to each mesh surface.
double d[n];
for (int k = 0; k < n; ++k) {
if (std::fabs(u[k]) < FP_PRECISION) {
@ -506,14 +540,13 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
}
}
// Assign the next tally bin and the score.
// Pick the closest mesh surface and append this traversal to the output.
auto j = std::min_element(d, d+n) - d;
double distance = d[j];
bins.push_back(get_bin_from_indices(ijk0));
lengths.push_back(distance / total_distance);
// Translate the starting coordintes by the distance to the oncoming mesh
// surface.
// Translate to the oncoming mesh surface.
r0 += distance * u;
// Increment the indices into the next mesh cell.
@ -536,7 +569,8 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
}
}
void RegularMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins) const
void RegularMesh::surface_bins_crossed(const Particle* p,
std::vector<int>& bins) const
{
// ========================================================================
// Determine if the track intersects the tally mesh.
@ -555,12 +589,15 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins
get_indices(r1, ijk1, &end_in_mesh);
// Check if the track intersects any part of the mesh.
if (!start_in_mesh && !end_in_mesh) {
if (!intersects(r0, r1)) return;
if (!start_in_mesh) {
Position r0_copy = r0;
int ijk0_copy[n];
for (int i = 0; i < n; ++i) ijk0_copy[i] = ijk0[i];
if (!intersects(r0_copy, r1, ijk0_copy)) return;
}
// ========================================================================
// Figure out which mesh cell to tally.
// Find which mesh surfaces are crossed.
// Calculate number of surface crossings
int n_cross = 0;
@ -587,7 +624,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins
double distance = INFTY;
for (int i = 0; i < n; ++i) {
if (u[i] == 0) {
d[i] = INFINITY;
d[i] = INFTY;
} else {
d[i] = (xyz_cross[i] - r0[i])/u[i];
}