crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Update get_bin interface

Browse source
This commit is contained in:
Joffrey Dorville 2026-06-01 16:29:48 -05:00
parent 5fdb700e73
commit 6eef5f48a7
10 changed files with 56 additions and 48 deletions
Download patch file Download diff file Expand all files Collapse all files

32
include/openmc/mesh.h
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,9 @@ 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 +648,9 @@ 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 +839,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 +951,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 +1006,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 +1096,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

2
include/openmc/random_ray/flat_source_domain.h
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;
//----------------------------------------------------------------------------

2
src/cmfd_solver.cpp
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) {

2
src/eigenvalue.cpp
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!");

48
src/mesh.cpp
View file

@ -1023,12 +1023,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;
@ -1066,11 +1066,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;
@ -1104,7 +1104,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) {
@ -1171,7 +1171,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
@ -1249,7 +1249,7 @@ 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);
@ -1586,7 +1586,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) {
@ -1827,7 +1827,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);
@ -1843,7 +1843,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]);
@ -2042,7 +2042,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;
}
@ -2120,7 +2120,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);
@ -2137,7 +2137,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]);
@ -2371,7 +2371,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;
}
@ -2656,13 +2656,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);
}
}
}
@ -3095,7 +3096,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);
@ -3114,7 +3115,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;
@ -3134,7 +3135,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);
@ -3227,7 +3228,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) {
@ -3311,9 +3312,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;
}
@ -3866,7 +3868,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);
@ -3967,7 +3969,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);

9
src/random_ray/flat_source_domain.cpp
View file

@ -1600,7 +1600,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;
@ -1825,18 +1825,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;
}

2
src/tallies/filter_mesh.cpp
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);

3
src/tallies/filter_meshborn.cpp
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);

2
src/tallies/filter_meshmaterial.cpp
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()) {

2
src/weight_windows.cpp
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)