This commit is contained in:
Joffrey Dorville 2026-07-19 23:23:14 -05:00 committed by GitHub
commit f183512eaa
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
15 changed files with 184 additions and 54 deletions

View file

@ -196,8 +196,9 @@ public:
//! Get bin at a given position in space
//
//! \param[in] r Position to get bin for
//! \param[in] u Direction of the particle
//! \return Mesh bin
virtual int get_bin(Position r) const = 0;
virtual int get_bin(Position r, Direction u) const = 0;
//! Get the number of mesh cells.
virtual int n_bins() const = 0;
@ -321,7 +322,7 @@ public:
virtual Position sample_element(const MeshIndex& ijk, uint64_t* seed) const;
int get_bin(Position r) const override;
int get_bin(Position r, Direction u) const override;
int n_bins() const override;
@ -360,9 +361,10 @@ public:
//! Get mesh indices given a position
//
//! \param[in] r Position to get indices for
//! \param[in] u Direction of the particle
//! \param[out] in_mesh Whether position is in mesh
//! \return Array of mesh indices
virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
virtual MeshIndex get_indices(Position r, Direction u, bool& in_mesh) const;
//! Get mesh indices corresponding to a mesh bin
//
@ -373,8 +375,9 @@ public:
//! Get mesh index in a particular direction
//!
//! \param[in] r Coordinate to get index for
//! \param[in] u Direction of the particle
//! \param[in] i Direction index
virtual int get_index_in_direction(double r, int i) const = 0;
virtual int get_index_in_direction(double r, double u, int i) const = 0;
//! Get the coordinate for the mesh grid boundary in the positive direction
//!
@ -484,7 +487,7 @@ public:
RegularMesh(hid_t group);
// Overridden methods
int get_index_in_direction(double r, int i) const override;
int get_index_in_direction(double r, double u, int i) const override;
virtual std::string get_mesh_type() const override;
@ -537,7 +540,7 @@ public:
RectilinearMesh(hid_t group);
// Overridden methods
int get_index_in_direction(double r, int i) const override;
int get_index_in_direction(double r, double u, int i) const override;
virtual std::string get_mesh_type() const override;
@ -580,9 +583,10 @@ public:
CylindricalMesh(hid_t group);
// Overridden methods
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
virtual MeshIndex get_indices(
Position r, Direction u, bool& in_mesh) const override;
int get_index_in_direction(double r, int i) const override;
int get_index_in_direction(double r, double u, int i) const override;
virtual std::string get_mesh_type() const override;
@ -645,9 +649,10 @@ public:
SphericalMesh(hid_t group);
// Overridden methods
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
virtual MeshIndex get_indices(
Position r, Direction u, bool& in_mesh) const override;
int get_index_in_direction(double r, int i) const override;
int get_index_in_direction(double r, double u, int i) const override;
virtual std::string get_mesh_type() const override;
@ -836,7 +841,7 @@ public:
void bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins, vector<double>& lengths) const override;
int get_bin(Position r) const override;
int get_bin(Position r, Direction u) const override;
int n_bins() const override;
@ -948,8 +953,9 @@ private:
//! Get the mesh cell index for a given position
//
//! \param[in] r Position to get index for
//! \param[in] u Direction of the particle
//! \param[in,out] in_mesh Whether position is in the mesh
int get_index(const Position& r, bool* in_mesh) const;
int get_index(const Position& r, const Direction& u, bool* in_mesh) const;
//! Get the mesh cell index from a bin
//
@ -1002,7 +1008,7 @@ public:
Position sample_element(int32_t bin, uint64_t* seed) const override;
virtual int get_bin(Position r) const override;
virtual int get_bin(Position r, Direction u) const override;
int n_bins() const override;
@ -1092,7 +1098,7 @@ public:
void write(const std::string& filename) const override;
int get_bin(Position r) const override;
int get_bin(Position r, Direction u) const override;
protected:
// Overridden methods

View file

@ -352,6 +352,10 @@ public:
Position& r_born() { return r_born_; }
const Position& r_born() const { return r_born_; }
// Direction at birth
Direction& u_born() { return u_born_; }
const Direction& u_born() const { return u_born_; }
// Coordinates of last collision or reflective/periodic surface
// crossing for current tallies
Position& r_last_current() { return r_last_current_; }
@ -426,6 +430,7 @@ private:
vector<int> cell_last_; //!< coordinates for all levels
Position r_born_; //!< coordinates at birth
Position u_born_; //!< direction at birth
Position r_last_current_; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies

View file

@ -70,7 +70,7 @@ public:
}
int64_t lookup_base_source_region_idx(const GeometryState& p) const;
SourceRegionKey lookup_source_region_key(const GeometryState& p) const;
int64_t lookup_mesh_bin(int64_t sr, Position r) const;
int64_t lookup_mesh_bin(int64_t sr, Position r, const Direction& u) const;
int lookup_mesh_idx(int64_t sr) const;
//----------------------------------------------------------------------------

View file

@ -94,7 +94,7 @@ tensor::Tensor<double> count_bank_sites(
const auto& site = simulation::source_bank[i];
// determine scoring bin for CMFD mesh
int mesh_bin = cmfd::mesh->get_bin(site.r);
int mesh_bin = cmfd::mesh->get_bin(site.r, site.u);
// if outside mesh, skip particle
if (mesh_bin < 0) {

View file

@ -660,7 +660,7 @@ void ufs_count_sites()
double ufs_get_weight(const Particle& p)
{
// Determine indices on ufs mesh for current location
int mesh_bin = simulation::ufs_mesh->get_bin(p.r());
int mesh_bin = simulation::ufs_mesh->get_bin(p.r(), p.u());
if (mesh_bin < 0) {
p.write_restart();
fatal_error("Source site outside UFS mesh!");

View file

@ -1024,12 +1024,12 @@ ElementType UnstructuredMesh::element_type(int bin) const
}
StructuredMesh::MeshIndex StructuredMesh::get_indices(
Position r, bool& in_mesh) const
Position r, Direction u, bool& in_mesh) const
{
MeshIndex ijk;
in_mesh = true;
for (int i = 0; i < n_dimension_; ++i) {
ijk[i] = get_index_in_direction(r[i], i);
ijk[i] = get_index_in_direction(r[i], u[i], i);
if (ijk[i] < 1 || ijk[i] > shape_[i])
in_mesh = false;
@ -1067,11 +1067,11 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const
return ijk;
}
int StructuredMesh::get_bin(Position r) const
int StructuredMesh::get_bin(Position r, Direction u) const
{
// Determine indices
bool in_mesh;
MeshIndex ijk = get_indices(r, in_mesh);
MeshIndex ijk = get_indices(r, u, in_mesh);
if (!in_mesh)
return -1;
@ -1118,7 +1118,7 @@ tensor::Tensor<double> StructuredMesh::count_sites(
const auto& site = bank[i];
// determine scoring bin for entropy mesh
int mesh_bin = get_bin(site.r);
int mesh_bin = get_bin(site.r, site.u);
// if outside mesh, skip particle
if (mesh_bin < 0) {
@ -1185,7 +1185,7 @@ void StructuredMesh::raytrace_mesh(
// Calculate index of current cell. Offset the position a tiny bit in
// direction of flight
MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh);
MeshIndex ijk = get_indices(global_r + TINY_BIT * u, u, in_mesh);
// if track is very short, assume that it is completely inside one cell.
// Only the current cell will score and no surfaces
@ -1263,7 +1263,8 @@ void StructuredMesh::raytrace_mesh(
// Calculate the new cell index and update all distances to next
// surfaces.
ijk = get_indices(global_r + (traveled_distance + TINY_BIT) * u, in_mesh);
ijk =
get_indices(global_r + (traveled_distance + TINY_BIT) * u, u, in_mesh);
for (int k = 0; k < n; ++k) {
distances[k] =
distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);
@ -1484,9 +1485,18 @@ RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group}
}
}
int RegularMesh::get_index_in_direction(double r, int i) const
int RegularMesh::get_index_in_direction(double r, double u, int i) const
{
return std::ceil((r - lower_left_[i]) / width_[i]);
int idx = std::ceil((r - lower_left_[i]) / width_[i]);
// If on upper boundary with positive direction, use next index
if (r == lower_left_[i] + width_[i] * idx) {
if (u > 0.0) {
idx++;
}
}
return idx;
}
const std::string RegularMesh::mesh_type = "regular";
@ -1591,7 +1601,7 @@ tensor::Tensor<double> RegularMesh::count_sites(
const auto& site = bank[i];
// determine scoring bin for entropy mesh
int mesh_bin = get_bin(site.r);
int mesh_bin = get_bin(site.r, site.u);
// if outside mesh, skip particle
if (mesh_bin < 0) {
@ -1725,9 +1735,25 @@ int RectilinearMesh::set_grid()
return 0;
}
int RectilinearMesh::get_index_in_direction(double r, int i) const
int RectilinearMesh::get_index_in_direction(double r, double u, int i) const
{
return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
int idx = lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
// If on lower boundary with negative direction, use previous index
if (r == grid_[i][idx - 1]) {
if (u < 0) {
idx--;
}
}
// If on upper boundary with positive direction, use next index
if (r == grid_[i][idx]) {
if (u > 0) {
idx++;
}
}
return idx;
}
std::pair<vector<double>, vector<double>> RectilinearMesh::plot(
@ -1816,7 +1842,7 @@ std::string CylindricalMesh::get_mesh_type() const
}
StructuredMesh::MeshIndex CylindricalMesh::get_indices(
Position r, bool& in_mesh) const
Position r, Direction u, bool& in_mesh) const
{
r = local_coords(r);
@ -1832,7 +1858,7 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices(
mapped_r[1] += 2 * M_PI;
}
MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
MeshIndex idx = StructuredMesh::get_indices(mapped_r, u, in_mesh);
idx[1] = sanitize_phi(idx[1]);
@ -2031,7 +2057,7 @@ int CylindricalMesh::set_grid()
return 0;
}
int CylindricalMesh::get_index_in_direction(double r, int i) const
int CylindricalMesh::get_index_in_direction(double r, double u, int i) const
{
return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
}
@ -2109,7 +2135,7 @@ std::string SphericalMesh::get_mesh_type() const
}
StructuredMesh::MeshIndex SphericalMesh::get_indices(
Position r, bool& in_mesh) const
Position r, Direction u, bool& in_mesh) const
{
r = local_coords(r);
@ -2126,7 +2152,7 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(
mapped_r[2] += 2 * M_PI;
}
MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
MeshIndex idx = StructuredMesh::get_indices(mapped_r, u, in_mesh);
idx[1] = sanitize_theta(idx[1]);
idx[2] = sanitize_phi(idx[2]);
@ -2360,7 +2386,7 @@ int SphericalMesh::set_grid()
return 0;
}
int SphericalMesh::get_index_in_direction(double r, int i) const
int SphericalMesh::get_index_in_direction(double r, double u, int i) const
{
return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
}
@ -2645,13 +2671,14 @@ extern "C" int openmc_mesh_get_plot_bins(int32_t index, Position origin,
#pragma omp parallel
{
Position r = xyz;
Direction placeholder_u = Direction(1.0, 0.0, 0.0);
#pragma omp for
for (int y = 0; y < pixel_height; y++) {
r[out_i] = xyz[out_i] - out_pixel * y;
for (int x = 0; x < pixel_width; x++) {
r[in_i] = xyz[in_i] + in_pixel * x;
data[pixel_width * y + x] = mesh->get_bin(r);
data[pixel_width * y + x] = mesh->get_bin(r, placeholder_u);
}
}
}
@ -3084,7 +3111,7 @@ void MOABMesh::bins_crossed(Position r0, Position r1, const Direction& u,
// score to that tet and return.
if (hits.size() == 0) {
Position midpoint = r0 + u * (track_len * 0.5);
int bin = this->get_bin(midpoint);
int bin = this->get_bin(midpoint, u);
if (bin != -1) {
bins.push_back(bin);
lengths.push_back(1.0);
@ -3103,7 +3130,7 @@ void MOABMesh::bins_crossed(Position r0, Position r1, const Direction& u,
// find the midpoint of this segment
Position midpoint = current + u * (segment_length * 0.5);
// try to find a tet for this position
int bin = this->get_bin(midpoint);
int bin = this->get_bin(midpoint, u);
// determine the start point for this segment
current = r0 + u * hit;
@ -3123,7 +3150,7 @@ void MOABMesh::bins_crossed(Position r0, Position r1, const Direction& u,
Position segment_start = r0 + u * hits.back();
double segment_length = track_len - hits.back();
Position midpoint = segment_start + u * (segment_length * 0.5);
int bin = this->get_bin(midpoint);
int bin = this->get_bin(midpoint, u);
if (bin != -1) {
bins.push_back(bin);
lengths.push_back(segment_length / track_len);
@ -3216,7 +3243,7 @@ double MOABMesh::tet_volume(moab::EntityHandle tet) const
return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0]));
}
int MOABMesh::get_bin(Position r) const
int MOABMesh::get_bin(Position r, Direction u) const
{
moab::EntityHandle tet = get_tet(r);
if (tet == 0) {
@ -3300,9 +3327,10 @@ int MOABMesh::get_bin_from_index(int idx) const
return ehs_[idx] - ehs_[0];
}
int MOABMesh::get_index(const Position& r, bool* in_mesh) const
int MOABMesh::get_index(
const Position& r, const Direction& u, bool* in_mesh) const
{
int bin = get_bin(r);
int bin = get_bin(r, u);
*in_mesh = bin != -1;
return bin;
}
@ -3855,7 +3883,7 @@ void LibMesh::bins_crossed(Position r0, Position r1, const Direction& u,
fatal_error("Tracklength tallies on libMesh instances are not implemented.");
}
int LibMesh::get_bin(Position r) const
int LibMesh::get_bin(Position r, Direction u) const
{
// look-up a tet using the point locator
libMesh::Point p(r.x, r.y, r.z);
@ -3956,7 +3984,7 @@ void AdaptiveLibMesh::write(const std::string& filename) const
this->id_));
}
int AdaptiveLibMesh::get_bin(Position r) const
int AdaptiveLibMesh::get_bin(Position r, Direction u) const
{
// look-up a tet using the point locator
libMesh::Point p(r.x, r.y, r.z);

View file

@ -165,6 +165,7 @@ void Particle::from_source(const SourceSite* src)
r() = src->r;
u() = src->u;
r_born() = src->r;
u_born() = src->u;
r_last_current() = src->r;
r_last() = src->r;
u_last() = src->u;

View file

@ -1608,7 +1608,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
}
} else {
Mesh* mesh = model::meshes[mesh_idx].get();
int bin_found = mesh->get_bin(r + TINY_BIT * u);
int bin_found = mesh->get_bin(r + TINY_BIT * u, u);
if (bin_found != sr_key.mesh_bin) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
@ -1833,18 +1833,19 @@ SourceRegionKey FlatSourceDomain::lookup_source_region_key(
const GeometryState& gs) const
{
int64_t sr = lookup_base_source_region_idx(gs);
int64_t mesh_bin = lookup_mesh_bin(sr, gs.r());
int64_t mesh_bin = lookup_mesh_bin(sr, gs.r(), gs.u());
return SourceRegionKey {sr, mesh_bin};
}
// Determines the mesh bin that corresponds to a particular base source region
// index and position.
int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const
int64_t FlatSourceDomain::lookup_mesh_bin(
int64_t sr, Position r, const Direction& u) const
{
int mesh_idx = lookup_mesh_idx(sr);
int mesh_bin = 0;
if (mesh_idx != C_NONE) {
mesh_bin = model::meshes[mesh_idx]->get_bin(r);
mesh_bin = model::meshes[mesh_idx]->get_bin(r, u);
}
return mesh_bin;
}

View file

@ -58,7 +58,7 @@ void MeshFilter::get_all_bins(
}
if (estimator != TallyEstimator::TRACKLENGTH) {
auto bin = model::meshes[mesh_]->get_bin(r);
auto bin = model::meshes[mesh_]->get_bin(r, u);
if (bin >= 0) {
match.bins_.push_back(bin);
match.weights_.push_back(1.0);

View file

@ -11,13 +11,14 @@ void MeshBornFilter::get_all_bins(
const Particle& p, TallyEstimator estimator, FilterMatch& match) const
{
Position r_born = p.r_born();
Direction u_born = p.u_born();
// apply translation if present
if (translated_) {
r_born -= translation();
}
auto bin = model::meshes[mesh_]->get_bin(r_born);
auto bin = model::meshes[mesh_]->get_bin(r_born, u_born);
if (bin >= 0) {
match.bins_.push_back(bin);
match.weights_.push_back(1.0);

View file

@ -122,7 +122,7 @@ void MeshMaterialFilter::get_all_bins(
}
if (estimator != TallyEstimator::TRACKLENGTH) {
int32_t index_element = model::meshes[mesh_]->get_bin(r);
int32_t index_element = model::meshes[mesh_]->get_bin(r, u);
if (index_element >= 0) {
auto search = map_.find({index_element, p.material()});
if (search != map_.end()) {

View file

@ -290,7 +290,7 @@ std::pair<bool, WeightWindow> WeightWindows::get_weight_window(
// Get mesh index for particle's position
const auto& mesh = this->mesh();
int mesh_bin = mesh->get_bin(p.r());
int mesh_bin = mesh->get_bin(p.r(), p.u());
// particle is outside the weight window mesh
if (mesh_bin < 0)

View file

@ -255,3 +255,91 @@ TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical")
REQUIRE(regular_mesh_hdf5->lower_left() == regular_mesh_xml->lower_left());
REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right());
}
TEST_CASE("Test get_index_in_direction - regular")
{
// The XML data as a string
std::string xml_string = R"(
<mesh id="1">
<dimension>2 1 1</dimension>
<lower_left>-2 -2 -2</lower_left>
<upper_right>2 2 2</upper_right>
</mesh>
)";
// Create a pugixml document object
pugi::xml_document doc;
// Load the XML from the string
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
pugi::xml_node root = doc.child("mesh");
auto mesh = RegularMesh(root);
REQUIRE(mesh.get_index_in_direction(-4.01, 1.0, 0) == -1);
REQUIRE(mesh.get_index_in_direction(-4.0, -1.0, 0) == -1);
REQUIRE(mesh.get_index_in_direction(-4.0, 1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-3.99, 1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-2.01, 1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-2.0, -1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-2.0, 1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(-1.99, 1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(-0.01, 1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(0.0, -1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(0.0, 1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(0.01, 1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(1.99, -1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(2.0, -1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(2.0, 1.0, 0) == 3);
REQUIRE(mesh.get_index_in_direction(2.01, -1.0, 0) == 3);
REQUIRE(mesh.get_index_in_direction(3.9, -1.0, 0) == 3);
REQUIRE(mesh.get_index_in_direction(4.0, -1.0, 0) == 3);
REQUIRE(mesh.get_index_in_direction(4.0, 1.0, 0) == 4);
REQUIRE(mesh.get_index_in_direction(4.1, -1.0, 0) == 4);
}
TEST_CASE("Test get_index_in_direction - rectilinear")
{
// The XML data as a string
std::string xml_string = R"(
<mesh id="1" type="rectilinear">
<x_grid>-1.0 0.0 2.0</x_grid>
<y_grid>-1.0 1.0</y_grid>
<z_grid>-1.0 1.0</z_grid>
</mesh>
)";
// Create a pugixml document object
pugi::xml_document doc;
// Load the XML from the string
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
pugi::xml_node root = doc.child("mesh");
auto mesh = RectilinearMesh(root);
REQUIRE(mesh.get_index_in_direction(-5.0, -1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-1.1, 1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-1.0, -1.0, 0) == 0);
REQUIRE(mesh.get_index_in_direction(-1.0, 1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(-0.9, -1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(-0.1, 1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(0.0, -1.0, 0) == 1);
REQUIRE(mesh.get_index_in_direction(0.0, 1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(0.1, 1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(1.99, -1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(2.0, -1.0, 0) == 2);
REQUIRE(mesh.get_index_in_direction(2.0, 1.0, 0) == 3);
REQUIRE(mesh.get_index_in_direction(2.1, -1.0, 0) == 3);
REQUIRE(mesh.get_index_in_direction(5.0, -1.0, 0) == 3);
}

View file

@ -1 +1 @@
386e507008ed3c72c6e1a101aafc01cfaaff2c0b10555558d1eab6332441f7b17264b55002d721151bac2f3e7c1a8aac4c5ed243f5270533d171850f985206ed
0d7b17d4e364bda2be21feb2e5b2aae4be69fc8ed634c4f45062e8efc61b7bd24dcf448837692fd603afe3756501fe8821d134f2d6289179618f85178ddb8793

View file

@ -1 +1 @@
d17a437262d3316985fba4b48e21a7fcd5f32ac2d96ef0e66849f567deaeadc1e0f18d5d0bf5e9cab4246dbf823e7f43f249a1f67e928b27ae8c70b89f3cefbf
2f1c8efbf6ab0b0c7319fadd7ee8f27b545c1ebab94ebe15cf1133b54227487542966b6d798970c8b37bc19e7041e687687d5285803241bb83f4c8dc40a4d276