diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0d8189caa1..e7af15cb08 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -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& bins, vector& 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 diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 6f51af34de..1b89f40254 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -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; //---------------------------------------------------------------------------- diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 714a5bf3ac..1b038e3ba8 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -94,7 +94,7 @@ tensor::Tensor 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) { diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index bc8dcc9ea7..dcae75801c 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -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!"); diff --git a/src/mesh.cpp b/src/mesh.cpp index 90b269a6ab..6f47d0973b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -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 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 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); diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 06c6ef14d7..c43d9a0d15 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -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; } diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index a0698992d0..c9e11072ef 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -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); diff --git a/src/tallies/filter_meshborn.cpp b/src/tallies/filter_meshborn.cpp index c95dc3dc78..24203bc60c 100644 --- a/src/tallies/filter_meshborn.cpp +++ b/src/tallies/filter_meshborn.cpp @@ -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); diff --git a/src/tallies/filter_meshmaterial.cpp b/src/tallies/filter_meshmaterial.cpp index b45bb4164e..5c31c271ac 100644 --- a/src/tallies/filter_meshmaterial.cpp +++ b/src/tallies/filter_meshmaterial.cpp @@ -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()) { diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 0614110cd3..9b59fd738a 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -290,7 +290,7 @@ std::pair 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)