From fa3fc916816fb79c19f95749aec12fdf32e41f53 Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Tue, 16 Jun 2026 16:05:43 -0400 Subject: [PATCH 1/7] angular flux init --- include/openmc/random_ray/source_region.h | 22 ++++++++++++++++++++++ src/random_ray/source_region.cpp | 1 + 2 files changed, 23 insertions(+) diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 65f2e6bc41c..f231b30fee2 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -180,6 +180,7 @@ class SourceRegionHandle { float* source_; float* external_source_; double* scalar_flux_final_; + vector* angular_flux_new_; MomentArray* source_gradients_; MomentArray* flux_moments_old_; @@ -279,6 +280,9 @@ class SourceRegionHandle { double& scalar_flux_final(int g) { return scalar_flux_final_[g]; } const double scalar_flux_final(int g) const { return scalar_flux_final_[g]; } + double& angular_flux_new(int g, int angle) { return angular_flux_new_[g]; } + const double angular_flux_new(int g) const { return angular_flux_new_[g]; } + float& source(int g) { return source_[g]; } const float source(int g) const { return source_[g]; } @@ -375,6 +379,9 @@ class SourceRegion { //!< active iterations (used for plotting, //!< or computing adjoint sources) + vector + angular_flux_new_; //!< The angular flux from the current iteration + vector source_gradients_; //!< The linear source gradients vector flux_moments_old_; //!< The linear flux moments from the previous iteration @@ -572,6 +579,20 @@ class SourceRegionContainer { return scalar_flux_final_[se]; } + double& angular_flux_new(int64_t sr, int g) + { + return angular_flux_new_[index(sr, g)]; + } + const double angular_flux_new(int64_t sr, int g) const + { + return angular_flux_new_[index(sr, g)]; + } + double& angular_flux_new(int64_t se) { return angular_flux_new_[se]; } + const double angular_flux_new(int64_t se) const + { + return angular_flux_new_[se]; + } + float& source(int64_t sr, int g) { return source_[index(sr, g)]; } const float source(int64_t sr, int g) const { return source_[index(sr, g)]; } float& source(int64_t se) { return source_[se]; } @@ -671,6 +692,7 @@ class SourceRegionContainer { vector scalar_flux_old_; vector scalar_flux_new_; vector scalar_flux_final_; + vector angular_flux_new_; vector source_; vector external_source_; diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 78543c5ab53..01d8d6b6962 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -26,6 +26,7 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) scalar_flux_new_(sr.scalar_flux_new_.data()), source_(sr.source_.data()), external_source_(sr.external_source_.data()), scalar_flux_final_(sr.scalar_flux_final_.data()), + angular_flux_new_(sr.angular_flux_new_.data()), source_gradients_(sr.source_gradients_.data()), flux_moments_old_(sr.flux_moments_old_.data()), flux_moments_new_(sr.flux_moments_new_.data()), From 8ff54d8175c572b3d201c5935cfa70e0e8384d76 Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Thu, 13 Aug 2026 12:06:17 -0400 Subject: [PATCH 2/7] Particle, RandomRay, SourceRegion and TallyTask changes to support angular flux accumulation --- include/openmc/particle.h | 4 ++ include/openmc/random_ray/random_ray.h | 3 + include/openmc/random_ray/source_region.h | 60 +++++++++++----- src/boundary_condition.cpp | 3 + src/random_ray/random_ray.cpp | 86 ++++++++++++++++++++--- src/random_ray/source_region.cpp | 28 ++++++-- 6 files changed, 151 insertions(+), 33 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 8db9721baa8..d4191213f2d 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -106,6 +106,10 @@ class Particle : public ParticleData { void cross_periodic_bc( const Surface& surf, Position new_r, Direction new_u, int new_surface); + //! Reset angular flux tally bin if in RandomRay mode + //! (does nothing for MC particles) + virtual void direction_changed() {} + //! mark a particle as lost and create a particle restart file //! \param message A warning message to display virtual void mark_as_lost(const char* message) override; diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index b61d2d67aa8..32a5657ed27 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -36,12 +36,14 @@ class RandomRay : public Particle { SourceRegionHandle& srh, double distance, bool is_active, Position r); void attenuate_flux_linear_source_void( SourceRegionHandle& srh, double distance, bool is_active, Position r); + void direction_changed() override { angular_bin_ = C_NONE; }; void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain); uint64_t transport_history_based_single_ray(); SourceSite sample_prng(); SourceSite sample_halton(); SourceSite sample_s2(); + int angular_bin(); // accessor for current angular quadrature bin index //---------------------------------------------------------------------------- // Static data members @@ -67,6 +69,7 @@ class RandomRay : public Particle { int negroups_; int ntemperature_; + int angular_bin_ {C_NONE}; FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source // data needed for ray transport double distance_travelled_ {0}; diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index f231b30fee2..db28adda0d6 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -54,9 +54,17 @@ struct TallyTask { int filter_idx; int score_idx; int score_type; - TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type) + // Angular quadrature bin index this task applies to. C_NONE (the default) + // means the task is angle-independent and scores against the source + // region's scalar flux, as normal. A non-negative value means the task + // is scoring the angular flux tally for that specific angular bin + // (source_regions_.angular_flux_new(sr, g, angle_bin)), and is only ever + // created for source regions with an external source present. + int angle_bin {C_NONE}; + TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type, + int angle_bin = C_NONE) : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx), - score_type(score_type) + score_type(score_type), angle_bin(angle_bin) {} TallyTask() = default; @@ -65,7 +73,8 @@ struct TallyTask { bool operator==(const TallyTask& other) const { return tally_idx == other.tally_idx && filter_idx == other.filter_idx && - score_idx == other.score_idx && score_type == other.score_type; + score_idx == other.score_idx && score_type == other.score_type && + angle_bin == other.angle_bin; } struct HashFunctor { @@ -76,6 +85,7 @@ struct TallyTask { hash_combine(seed, task.filter_idx); hash_combine(seed, task.score_idx); hash_combine(seed, task.score_type); + hash_combine(seed, task.angle_bin); return seed; } }; @@ -141,6 +151,7 @@ class SourceRegionHandle { //---------------------------------------------------------------------------- // Public Data members int negroups_; + int nangles_ {1}; //!< Number of angular bins for angular flux binning bool is_numerical_fp_artifact_ {false}; bool is_linear_ {false}; @@ -180,7 +191,7 @@ class SourceRegionHandle { float* source_; float* external_source_; double* scalar_flux_final_; - vector* angular_flux_new_; + double* angular_flux_new_; //!< only accumulated for ext. source biasing MomentArray* source_gradients_; MomentArray* flux_moments_old_; @@ -280,8 +291,14 @@ class SourceRegionHandle { double& scalar_flux_final(int g) { return scalar_flux_final_[g]; } const double scalar_flux_final(int g) const { return scalar_flux_final_[g]; } - double& angular_flux_new(int g, int angle) { return angular_flux_new_[g]; } - const double angular_flux_new(int g) const { return angular_flux_new_[g]; } + double& angular_flux_new(int g, int a) + { + return angular_flux_new_[g * nangles_ + a]; + } + const double angular_flux_new(int g, int a) const + { + return angular_flux_new_[g * nangles_ + a]; + } float& source(int g) { return source_[g]; } const float source(int g) const { return source_[g]; } @@ -319,7 +336,7 @@ class SourceRegion { public: //---------------------------------------------------------------------------- // Constructors - SourceRegion(int negroups, bool is_linear); + SourceRegion(int negroups, bool is_linear, int nangles = 1); SourceRegion() = default; //---------------------------------------------------------------------------- @@ -402,8 +419,8 @@ class SourceRegionContainer { public: //---------------------------------------------------------------------------- // Constructors - SourceRegionContainer(int negroups, bool is_linear) - : negroups_(negroups), is_linear_(is_linear) + SourceRegionContainer(int negroups, bool is_linear, int nangles = 1) + : negroups_(negroups), is_linear_(is_linear), nangles_(nangles) {} SourceRegionContainer() = default; @@ -579,18 +596,18 @@ class SourceRegionContainer { return scalar_flux_final_[se]; } - double& angular_flux_new(int64_t sr, int g) + double& angular_flux_new(int64_t sr, int g, int a) { - return angular_flux_new_[index(sr, g)]; + return angular_flux_new_[index_angle(sr, g, a)]; } - const double angular_flux_new(int64_t sr, int g) const + const double angular_flux_new(int64_t sr, int g, int a) const { - return angular_flux_new_[index(sr, g)]; + return angular_flux_new_[index_angle(sr, g, a)]; } - double& angular_flux_new(int64_t se) { return angular_flux_new_[se]; } - const double angular_flux_new(int64_t se) const + double& angular_flux_new(int64_t sea) { return angular_flux_new_[sea]; } + const double angular_flux_new(int64_t sea) const { - return angular_flux_new_[se]; + return angular_flux_new_[sea]; } float& source(int64_t sr, int g) { return source_[index(sr, g)]; } @@ -647,8 +664,14 @@ class SourceRegionContainer { void flux_swap(); int64_t n_source_regions() const { return n_source_regions_; } int64_t n_source_elements() const { return n_source_regions_ * negroups_; } + int64_t n_source_angular_elements() const + { + return n_source_regions_ * negroups_ * nangles_; + } int& negroups() { return negroups_; } const int negroups() const { return negroups_; } + int& nangles() { return nangles_; } + const int nangles() const { return nangles_; } bool& is_linear() { return is_linear_; } const bool is_linear() const { return is_linear_; } SourceRegionHandle get_source_region_handle(int64_t sr); @@ -659,6 +682,7 @@ class SourceRegionContainer { // Private Data Members int64_t n_source_regions_ {0}; int negroups_ {0}; + int nangles_ {1}; bool is_linear_ {false}; // SoA storage for scalar fields (one item per source region) @@ -713,6 +737,10 @@ class SourceRegionContainer { // Helper function for indexing inline int index(int64_t sr, int g) const { return sr * negroups_ + g; } + inline int64_t index_angle(int64_t sr, int g, int a) const + { + return (sr * negroups_ + g) * nangles_ + a; + } }; } // namespace openmc diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 5bbda483059..85adb200ee8 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -39,6 +39,7 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const { Direction u = surf.reflect(p.r(), p.u(), &p); u /= u.norm(); + p.direction_changed(); // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); @@ -54,6 +55,7 @@ void WhiteBC::handle_particle(Particle& p, const Surface& surf) const { Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed()); u /= u.norm(); + p.direction_changed(); // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); @@ -240,6 +242,7 @@ void RotationalPeriodicBC::handle_particle( new_u[zero_axis_idx_] = u[zero_axis_idx_]; new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_]; new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_]; + p.direction_changed(); // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index dde5023e44f..f88768984f2 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -457,8 +457,22 @@ void RandomRay::attenuate_flux_flat_source( if (is_active) { // Accumulate delta psi into new estimate of source region flux for // this iteration - for (int g = 0; g < negroups_; g++) { - srh.scalar_flux_new(g) += delta_psi_[g]; + if (srh.external_source_present()) { + // Accumulate angular flux data in regions with external source + // for later source biasing. + // + // TODO: replace the above condition with a different + // SourceRegionHandle flag to enable angular flux tallying more + // generally. + int a = angular_bin(); + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += delta_psi_[g]; + srh.angular_flux_new(g, a) += delta_psi_[g]; + } + } else { + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += delta_psi_[g]; + } } // Accomulate volume (ray distance) into this iteration's estimate @@ -498,8 +512,22 @@ void RandomRay::attenuate_flux_flat_source_void( // Accumulate delta psi into new estimate of source region flux for // this iteration - for (int g = 0; g < negroups_; g++) { - srh.scalar_flux_new(g) += angular_flux_[g] * distance; + if (srh.external_source_present()) { + // Accumulate angular flux data in regions with external source + // for later source biasing. + // + // TODO: replace the above condition with a different + // SourceRegionHandle flag to enable angular flux tallying more + // generally. + int a = angular_bin(); + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += angular_flux_[g] * distance; + srh.angular_flux_new(g, a) += angular_flux_[g] * distance; + } + } else { + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += angular_flux_[g] * distance; + } } // Accomulate volume (ray distance) into this iteration's estimate @@ -628,9 +656,24 @@ void RandomRay::attenuate_flux_linear_source( if (is_active) { // Accumulate deltas into the new estimate of source region flux for this // iteration - for (int g = 0; g < negroups_; g++) { - srh.scalar_flux_new(g) += delta_psi_[g]; - srh.flux_moments_new(g) += delta_moments_[g]; + if (srh.external_source_present()) { + // Accumulate angular flux data in regions with external source + // for later source biasing. + // + // TODO: replace the above condition with a different + // SourceRegionHandle flag to enable angular flux tallying more + // generally. + int a = angular_bin(); + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += delta_psi_[g]; + srh.flux_moments_new(g) += delta_moments_[g]; + srh.angular_flux_new(g, a) += delta_psi_[g]; + } + } else { + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += delta_psi_[g]; + srh.flux_moments_new(g) += delta_moments_[g]; + } } // Accumulate the volume (ray segment distance), centroid, and spatial @@ -732,9 +775,24 @@ void RandomRay::attenuate_flux_linear_source_void( // Accumulate delta psi into new estimate of source region flux for // this iteration, and update flux momements - for (int g = 0; g < negroups_; g++) { - srh.scalar_flux_new(g) += angular_flux_[g] * distance; - srh.flux_moments_new(g) += delta_moments_[g]; + if (srh.external_source_present()) { + // Accumulate angular flux data in regions with external source + // for later source biasing. + // + // TODO: replace the above condition with a different + // SourceRegionHandle flag to enable angular flux tallying more + // generally. + int a = angular_bin(); + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += angular_flux_[g] * distance; + srh.flux_moments_new(g) += delta_moments_[g]; + srh.angular_flux_new(g, a) += angular_flux_[g] * distance; + } + } else { + for (int g = 0; g < negroups_; g++) { + srh.scalar_flux_new(g) += angular_flux_[g] * distance; + srh.flux_moments_new(g) += delta_moments_[g]; + } } // Accumulate the volume (ray segment distance), centroid, and spatial @@ -900,4 +958,12 @@ SourceSite RandomRay::sample_s2() return site; } +int RandomRay::angular_bin() +{ + if (angular_bin_ == C_NONE) { + angular_bin_ = domain_->get_angular_bin(u()); + } + return angular_bin_; +} + } // namespace openmc diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 01d8d6b6962..5a305cfe27e 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -10,12 +10,16 @@ namespace openmc { // SourceRegionHandle implementation //============================================================================== SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) - : negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_), - temperature_idx_(&sr.temperature_idx_), density_mult_(&sr.density_mult_), - is_small_(&sr.is_small_), n_hits_(&sr.n_hits_), - is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_), - volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_), - volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_), + : negroups_(sr.scalar_flux_old_.size()), + nangles_(sr.angular_flux_new_.empty() + ? 1 + : sr.angular_flux_new_.size() / sr.scalar_flux_old_.size()), + material_(&sr.material_), temperature_idx_(&sr.temperature_idx_), + density_mult_(&sr.density_mult_), is_small_(&sr.is_small_), + n_hits_(&sr.n_hits_), is_linear_(sr.source_gradients_.size() > 0), + lock_(&sr.lock_), volume_(&sr.volume_), volume_t_(&sr.volume_t_), + volume_sq_(&sr.volume_sq_), volume_sq_t_(&sr.volume_sq_t_), + volume_naive_(&sr.volume_naive_), position_recorded_(&sr.position_recorded_), external_source_present_(&sr.external_source_present_), position_(&sr.position_), centroid_(&sr.centroid_), @@ -37,7 +41,7 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) //============================================================================== // SourceRegion implementation //============================================================================== -SourceRegion::SourceRegion(int negroups, bool is_linear) +SourceRegion::SourceRegion(int negroups, bool is_linear, int nangles) { if (settings::run_mode == RunMode::EIGENVALUE) { // If in eigenvalue mode, set starting flux to guess of 1 @@ -52,6 +56,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) scalar_flux_new_.assign(negroups, 0.0); source_.assign(negroups, 0.0); scalar_flux_final_.assign(negroups, 0.0); + angular_flux_new_.assign(negroups * nangles, 0.0f); tally_task_.resize(negroups); if (is_linear) { @@ -119,6 +124,11 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) // Tally tasks tally_task_.emplace_back(sr.tally_task_[g]); } + + // Angle- and energy-dependent flux + for (int ga = 0; ga < negroups_ * nangles_; ++ga) { + angular_flux_new_.push_back(sr.angular_flux_new_[ga]); + } } void SourceRegionContainer::assign( @@ -154,6 +164,7 @@ void SourceRegionContainer::assign( scalar_flux_old_.clear(); scalar_flux_new_.clear(); scalar_flux_final_.clear(); + angular_flux_new_.clear(); source_.clear(); external_source_.clear(); @@ -185,6 +196,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) { SourceRegionHandle handle; handle.negroups_ = negroups(); + handle.nangles_ = nangles(); handle.material_ = &material(sr); handle.temperature_idx_ = &temperature_idx(sr); handle.density_mult_ = &density_mult(sr); @@ -212,6 +224,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) handle.external_source_ = nullptr; } handle.scalar_flux_final_ = &scalar_flux_final(sr, 0); + handle.angular_flux_new_ = &angular_flux_new(sr, 0, 0); handle.tally_task_ = &tally_task(sr, 0); if (handle.is_linear_) { @@ -254,6 +267,7 @@ void SourceRegionContainer::adjoint_reset() std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0); } std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0); + std::fill(angular_flux_new_.begin(), angular_flux_new_.end(), 0.0f); std::fill(source_.begin(), source_.end(), 0.0f); std::fill(external_source_.begin(), external_source_.end(), 0.0f); std::fill(source_gradients_.begin(), source_gradients_.end(), From eae7736d5263f4c9bc6e642cf3e27ca0d9cf18dc Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Fri, 21 Aug 2026 14:03:16 -0400 Subject: [PATCH 3/7] angular flux tallying workflow init --- include/openmc/mesh.h | 77 +++++++ .../openmc/random_ray/flat_source_domain.h | 14 +- .../openmc/random_ray/linear_source_domain.h | 2 +- include/openmc/random_ray/random_ray.h | 2 +- include/openmc/tallies/filter.h | 1 + include/openmc/tallies/filter_meshangular.h | 41 ++++ openmc/filter.py | 10 +- openmc/lib/filter.py | 42 +++- openmc/lib/mesh.py | 5 +- openmc/mesh.py | 196 ++++++++++++++++++ src/mesh.cpp | 80 +++++++ src/random_ray/flat_source_domain.cpp | 44 +++- src/random_ray/linear_source_domain.cpp | 2 +- src/random_ray/random_ray.cpp | 2 +- src/random_ray/random_ray_simulation.cpp | 2 +- src/tallies/filter.cpp | 3 + src/tallies/filter_mesh.cpp | 6 +- src/tallies/filter_meshangular.cpp | 49 +++++ 18 files changed, 560 insertions(+), 18 deletions(-) create mode 100644 include/openmc/tallies/filter_meshangular.h create mode 100644 src/tallies/filter_meshangular.cpp diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3d11d39017a..9c10356b1d6 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -826,6 +826,83 @@ class UnstructuredMesh : public Mesh { virtual void initialize() = 0; }; +// Abstract class for meshes over directions (points on the unit sphere), +// as opposed to meshes over volumes. +class AngularMesh : public Mesh { +public: + AngularMesh() { n_dimension_ = 2; } + AngularMesh(pugi::xml_node node) : Mesh(node) { n_dimension_ = 2; } + AngularMesh(hid_t group) : Mesh(group) { n_dimension_ = 2; } + + // Angular meshes partition direction space, not a volume. + // Therefore, "bin crossing" methods are not used; only the "get_bin" method. + void bins_crossed(Position r0, Position r1, const Direction& u, + vector& bins, vector& lengths) const override + { + fatal_error("Angular meshes do not support spatial tracklength tallies."); + } + + void surface_bins_crossed(Position r0, Position r1, const Direction& u, + vector& bins) const override + { + fatal_error("Angular meshes do not support surface-crossing tallies."); + } + + int n_surface_bins() const override { return 0; } + + std::pair, vector> plot( + Position plot_ll, Position plot_ur) const override + { + return {{}, {}}; + } + + std::string bin_label(int bin) const override + { + return fmt::format("Element Index ({})", bin); + } + + Position lower_left() const override { return {-1., -1., -1.}; } + Position upper_right() const override { return {1., 1., 1.}; } +}; + +class UnitSpherePointset : public AngularMesh { +public: + UnitSpherePointset() = default; + explicit UnitSpherePointset(vector points); + UnitSpherePointset(pugi::xml_node node); + UnitSpherePointset(hid_t group); + + //! TODO: add sampling from within a spherical Voronoi cell + Position sample_element(int32_t bin, uint64_t* seed) const override + { + fatal_error( + "Sampling over a UnitSpherePointset angular mesh is not supported"); + } + + int get_bin(Direction u) const override; + + int n_bins() const override { return static_cast(points_.size()); } + + double volume(int bin) const override + { + fatal_error("Volume calculation over UnitSpherePointset is not supported"); + } + + void material_volumes(int nx, int ny, int nz, int max_materials, + int32_t* materials, double* volumes, double* bboxes) const override + { + fatal_error("material_volumes() is not supported for UnitSpherePointset") + } + + std::string get_mesh_type() const override { return mesh_type; } + static const std::string mesh_type; + + void to_hdf5_inner(hid_t group) const override; + + vector points_; + vector data_; +}; + #ifdef OPENMC_DAGMC_ENABLED class MOABMesh : public UnstructuredMesh { diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 09414fd4465..216fe0ec436 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -30,7 +30,7 @@ class FlatSourceDomain { virtual void update_single_neutron_source(SourceRegionHandle& srh); virtual void update_all_neutron_sources(); void compute_k_eff(); - virtual void normalize_scalar_flux_and_volumes( + virtual void normalize_flux_and_volumes( double total_active_distance_per_iteration); int64_t add_source_to_scalar_flux(); @@ -68,10 +68,16 @@ class FlatSourceDomain { { return source_regions_.n_source_regions() * negroups_; } + int64_t n_source_angular_elements() const + { + return source_regions_.n_source_angular_elements(); + } 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; int lookup_mesh_idx(int64_t sr) const; + int lookup_angular_bin(Direction u) const; + Direction angular_mesh_direction(int a) { return angular_bin_angles_[a]; } //---------------------------------------------------------------------------- // Static Data members @@ -174,6 +180,12 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Private data members int negroups_; // Number of energy groups in simulation + int nangles_; // Number of bins for any angular flux tallies + + vector angular_bin_angles_; // Directions corresponding to each + // bin midpoint in the angular flux + // tallying scheme, flattened to 1D: + // [x0, y0, z0, x1, y1, z1,...] double simulation_volume_; // Total physical volume of the simulation domain, as diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h index 0098c782001..2c083dcec0d 100644 --- a/include/openmc/random_ray/linear_source_domain.h +++ b/include/openmc/random_ray/linear_source_domain.h @@ -21,7 +21,7 @@ class LinearSourceDomain : public FlatSourceDomain { //---------------------------------------------------------------------------- // Methods void update_single_neutron_source(SourceRegionHandle& srh) override; - void normalize_scalar_flux_and_volumes( + void normalize_flux_and_volumes( double total_active_distance_per_iteration) override; void batch_reset() override; diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index 32a5657ed27..e4d0bb0d443 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -36,7 +36,7 @@ class RandomRay : public Particle { SourceRegionHandle& srh, double distance, bool is_active, Position r); void attenuate_flux_linear_source_void( SourceRegionHandle& srh, double distance, bool is_active, Position r); - void direction_changed() override { angular_bin_ = C_NONE; }; + void direction_changed() override { angular_bin_ = C_NONE; } void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain); uint64_t transport_history_based_single_ray(); diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 77b0d9f420d..d78ddadcf72 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -32,6 +32,7 @@ enum class FilterType { MATERIAL, MATERIALFROM, MESH, + MESH_ANGULAR, MESHBORN, MESH_MATERIAL, MESH_SURFACE, diff --git a/include/openmc/tallies/filter_meshangular.h b/include/openmc/tallies/filter_meshangular.h new file mode 100644 index 00000000000..8535f5bf699 --- /dev/null +++ b/include/openmc/tallies/filter_meshangular.h @@ -0,0 +1,41 @@ +#ifndef OPENMC_TALLIES_FILTER_MESHSURFACE_H +#define OPENMC_TALLIES_FILTER_MESHSURFACE_H + +#include "openmc/tallies/filter_mesh.h" + +namespace openmc { + +//============================================================================== +//! Indexes the direction of particle events to a mesh. +//============================================================================== + +class MeshAngularFilter : public MeshFilter { +public: + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "meshangular"; } + FilterType type() const override { return FilterType::MESH_ANGULAR; } + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; + + void to_statepoint(hid_t filter_group) const override; + + + //---------------------------------------------------------------------------- + // Accessors + + void set_translation(const Position& translation) const override + { + fatal_error("Angular mesh filters do not permit translation."); + } + + void set_translation(const double translation[3]) const override + { + fatal_error("Angular mesh filters do not permit translation."); + } +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_MESHSURFACE_H diff --git a/openmc/filter.py b/openmc/filter.py index 31d72b0dd4c..00a2f31014c 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -27,7 +27,7 @@ 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide', - 'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction', + 'weight', 'meshangular', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction', ) def _mesh_current_names(mesh): @@ -1349,6 +1349,14 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs): # Initialize a Pandas DataFrame from the mesh dictionary return pd.concat([df, pd.DataFrame(filter_dict)]) +class MeshAngularFilter(Filter): + """Bins tally events based on incident particle's direction, using + an angular mesh. + + """ + def __init__(self): + pass + class CollisionFilter(Filter): """Bins tally events based on the number of collisions. diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index cd011dc5d42..edda4ed77c6 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -20,12 +20,12 @@ 'Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter', 'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', - 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter', - 'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', - 'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter', - 'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', - 'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter', - 'ZernikeRadialFilter', 'filters' + 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshAngularFilter', + 'MeshBornFilter', 'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', + 'MuSurfaceFilter', 'ParentNuclideFilter', 'ParticleFilter', + 'ParticleProductionFilter', 'PolarFilter', 'ReactionFilter', + 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', + 'UniverseFilter', 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' ] # Tally functions @@ -106,6 +106,12 @@ c_int32, POINTER(c_double), c_size_t] _dll.openmc_mesh_filter_set_rotation.restype = c_int _dll.openmc_mesh_filter_set_rotation.errcheck = _error_handler +_dll.openmc_meshangular_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_meshangular_filter_get_mesh.restype = c_int +_dll.openmc_meshangular_filter_get_mesh.errcheck = _error_handler +_dll.openmc_meshangular_filter_set_mesh.argtypes = [c_int32, c_int32] +_dll.openmc_meshangular_filter_set_mesh.restype = c_int +_dll.openmc_meshangular_filter_set_mesh.errcheck = _error_handler _dll.openmc_meshborn_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_meshborn_filter_get_mesh.restype = c_int _dll.openmc_meshborn_filter_get_mesh.errcheck = _error_handler @@ -466,6 +472,29 @@ def rotation(self, rotation_data): _dll.openmc_mesh_filter_set_rotation( self._index, flat_rotation.ctypes.data_as(POINTER(c_double)), c_size_t(len(flat_rotation))) + +class MeshAngularFilter(MeshFilter): + """Angular mesh filter stored internally. + + """ + filter_type = 'meshangular' + + def __init__(self, mesh=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if mesh is not None: + self.mesh = mesh + + @mesh.setter + def mesh(self, mesh): + _dll.openmc_meshangular_filter_set_mesh(self._index, mesh._index) + + @property + def translation(self): + raise AttributeError("Angular mesh filters do not permit translation.") + + @translation.setter + def translation(self, translation): + raise AttributeError("Angular mesh filters do not permit translation.") class MeshBornFilter(Filter): """MeshBorn filter stored internally. @@ -725,6 +754,7 @@ class ZernikeRadialFilter(ZernikeFilter): 'material': MaterialFilter, 'materialfrom': MaterialFromFilter, 'mesh': MeshFilter, + 'meshangular': MeshAngularFilter, 'meshborn': MeshBornFilter, 'meshmaterial': MeshMaterialFilter, 'meshsurface': MeshSurfaceFilter, diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 19e6f74d7ad..f176422fbca 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -730,13 +730,16 @@ def set_grid(self, r_grid, theta_grid, phi_grid): class UnstructuredMesh(Mesh): pass +class UnitSpherePointset(Mesh): + pass _MESH_TYPE_MAP = { 'regular': RegularMesh, 'rectilinear': RectilinearMesh, 'cylindrical': CylindricalMesh, 'spherical': SphericalMesh, - 'unstructured': UnstructuredMesh + 'unstructured': UnstructuredMesh, + 'unitsphere_pointset': UnitSpherePointset } diff --git a/openmc/mesh.py b/openmc/mesh.py index 670fcceab67..9d73a6a4a87 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -307,6 +307,8 @@ def from_hdf5(cls, group: h5py.Group): return SphericalMesh.from_hdf5(group, mesh_id, mesh_name) elif mesh_type == 'unstructured': return UnstructuredMesh.from_hdf5(group, mesh_id, mesh_name) + elif mesh_type == 'angular_pointset': + return UnitSpherePointset.from_hdf5(group, mesh_id, mesh_name) else: raise ValueError('Unrecognized mesh type: "' + mesh_type + '"') @@ -354,6 +356,8 @@ def from_xml_element(cls, elem: ET.Element): mesh = SphericalMesh.from_xml_element(elem) elif mesh_type == 'unstructured': mesh = UnstructuredMesh.from_xml_element(elem) + elif mesh_type == 'angular_pointset': + mesh = UnitSpherePointset.from_xml_element(elem) else: raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.') @@ -3352,6 +3356,198 @@ def from_xml_element(cls, elem: ET.Element): return cls(filename, library, mesh_id, '', length_multiplier, options) +class AngularMesh(MeshBase): + """Base class for angular meshes of the unit sphere.""" + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + @property + @abstractmethod + def dimension(self): + pass + + @property + @abstractmethod + def n_dimension(self): + pass + + @property + @abstractmethod + def axis_labels(self): + """tuple of str : Names of the mesh axes, one per dimension.""" + pass + + @property + def n_elements(self): + pass + + # override non-applicable methods of MeshBase + def get_homogenized_materials(self, *args, **kwargs): + raise NotImplementedError( + "Material attributes are not available for angular meshes.") + + def material_volumes(self, *args, **kwargs): + raise NotImplementedError( + "Material attributes are not available for angular meshes.") + +class UnitSpherePointset(AngularMesh): + """Set of points on the unit sphere. Used for mesh-based angular flux + tallying based on a Voronoi diagram generated from the pointset. + + Parameters + ---------- + points : iterable of float + Unit-vector endpoints constituting the mesh. Must be either a numpy + array or a nested list and have shape (N,3), where each row provides + the [x, y, z] components of one such vector. + data : iterable of float + Scalar data (e.g. flux) assigned to each point in the mesh + mesh_id : int + Unique identifier for the mesh + name : str + Name of the mesh + + Attributes + ---------- + id : int + Unique identifier for the mesh + name : str + Name of the mesh + points : numpy array of float + Array storing the + data : numpy array of float + Scalar data assigned to each point in the mesh + """ + def __init__( + self, + points, + data: float | None = None, + mesh_id: int | None = None, + name: str = '', + ): + super().__init__(mesh_id, name) + self.points = points + self.data = data + + @property + def points(self): + return self._points + + @points.setter + def points(self, pts): + cv.check_type("unit vector pointset", pts, Iterable, Real) + pts = np.asarray(pts) + if pts.shape != (len(pts),3): + raise ValueError( + "Unit vector array for UnitSpherePointset must have shape (N,3).") + self._points = pts + + @property + def data(self): + return self._data + + @data.setter + def data(self, d): + if d is not None: + cv.check_type("AngularMesh data", d, Iterable, Real) + d = np.asarray(d).flatten() + if len(d) != len(self.points): + raise ValueError( + "Data for UnitSpherePointset does not match number of mesh elements.") + self._data = d + + @property + def n_elements(self): + return len(self.points) + + @property + def dimension(self): + return (self.n_elements,) + + @property + def n_dimension(self): + return 2 + + @property + def lower_left(self): + return np.array((-1., -1., -1.)) + + @property + def upper_right(self): + return np.array((1., 1., 1.)) + + @property + def axis_labels(self): + return ('element_index',) + + @property + def indices(self): + return [(i,) for i in range(self.n_elements)] + + @classmethod + def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): + points = np.asarray(group['points'][()]) + n = points.size // 3 + points = points.reshape(n, 3) + + data = group['data'][()] if 'data' in group else None + + return cls(points, data=data, mesh_id=mesh_id, name=name) + + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : lxml.etree._Element + XML element containing mesh data + + """ + element = super().to_xml_element() + element.set("type", "angular_pointset") + + # flatten to a (3*N,) array + pts = self.points.flatten() + + subelement = ET.SubElement(element, "points") + subelement.text = ' '.join(map(str, pts)) + + if self.data is not None: + subelement = ET.SubElement(element, "data") + subelement.text = ' '.join(map(str, self.data)) + + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element): + """Generate a unit sphere pointset from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.UnitSpherePointset + Unit-sphere pointset object + + """ + mesh_id = int(get_text(elem, 'id')) + + points = np.array(get_elem_list(elem, "points", float)) + n = points.size // 3 + points = points.reshape(n, 3) + + data_elem = elem.find("data") + if data_elem is not None: + data = get_elem_list(elem, "data", float) + else: + data = None + + return cls(points, data, mesh_id=mesh_id) + def _read_meshes(elem): """Generate dictionary of meshes from a given XML node diff --git a/src/mesh.cpp b/src/mesh.cpp index a0e497613c0..e288b455967 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2429,6 +2429,86 @@ double SphericalMesh::volume(const MeshIndex& ijk) const (std::cos(theta_i) - std::cos(theta_o)) * (phi_o - phi_i); } +//============================================================================== +// Angular mesh implementations +//============================================================================== + +const std::string UnitSpherePointset::mesh_type = "angular_pointset"; + +UnitSpherePointset::UnitSpherePointset(vector points) + : points_(std::move(points)) +{} + +UnitSpherePointset::UnitSpherePointset(pugi::xml_node node) : AngularMesh(node) +{ + if (check_for_node(node, "type")) { + auto temp = get_node_value(node, "type", true, true); + if (temp != mesh_type) + fatal_error(fmt::format("Invalid mesh type: {}", temp)); + } + + vector flat = get_node_array(node, "points"); + if (flat.size() % 3 != 0) { + fatal_error(fmt::format("Point array for unit sphere pointset mesh {} " + "does not contain a whole number of points.", + id_)); + } + int n = static_cast(flat.size() / 3); + points_.reserve(n); + for (int i = 0; i < n; ++i) + points_.push_back({flat[3 * i], flat[3 * i + 1], flat[3 * i + 2]}); +} + +UnitSpherePointset::UnitSpherePointset(hid_t group) : AngularMesh(group) +{ + if (object_exists(group, "type")) { + std::string temp; + read_dataset(group, "type", temp); + if (temp != mesh_type) + fatal_error(fmt::format("Invalid mesh type: {}", temp)); + } + + vector flat; + read_dataset(group, "points", flat); + int n = static_cast(flat.size() / 3); + points_.reserve(n); + for (int i = 0; i < n; ++i) + points_.push_back({flat[3 * i], flat[3 * i + 1], flat[3 * i + 2]}); + + if (object_exists(group, "data")) + read_dataset(group, "data", data_); +} + +void UnitSpherePointset::to_hdf5_inner(hid_t mesh_group) const +{ + int n = this->n_bins(); + vector flat(3 * n); + for (int i = 0; i < n; ++i) { + flat[3 * i + 0] = points_[i].x; + flat[3 * i + 1] = points_[i].y; + flat[3 * i + 2] = points_[i].z; + } + write_dataset(mesh_group, "points", flat); + + if (!data_.empty()) + write_dataset(mesh_group, "data", data_); +} + +int UnitSpherePointset::get_bin(Direction u) const +{ + int best = -1; + double best_dot = -2.0; + for (int i = 0; i < this->n_bins(); ++i) { + double d = points_[i].dot(u); + if (d > best_dot) { + best_dot = d; + best = i; + } + } + return best; +} + + //============================================================================== // Helper functions for the C API //============================================================================== diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 1346a2b23af..94376c9fefb 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -54,9 +54,13 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) } } + // Count the number of angular bins if performing angular flux tallying + // and store the angle set. + nangles = + // Initialize source regions. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; - source_regions_ = SourceRegionContainer(negroups_, is_linear); + source_regions_ = SourceRegionContainer(negroups_, is_linear, nangles_); // Initialize tally volumes if (volume_normalized_flux_tallies_) { @@ -92,6 +96,11 @@ void FlatSourceDomain::batch_reset() for (int64_t se = 0; se < n_source_elements(); se++) { source_regions_.scalar_flux_new(se) = 0.0; } + +#pragma omp parallel for + for (int64_t sea = 0; sea < n_source_angular_elements(); sea++) { + source_regions_.angular_flux_new(sea) = 0.0; + } } void FlatSourceDomain::accumulate_iteration_flux() @@ -165,7 +174,7 @@ void FlatSourceDomain::update_all_neutron_sources() } // Normalizes flux and updates simulation-averaged volume estimate -void FlatSourceDomain::normalize_scalar_flux_and_volumes( +void FlatSourceDomain::normalize_flux_and_volumes( double total_active_distance_per_iteration) { double normalization_factor = 1.0 / total_active_distance_per_iteration; @@ -179,6 +188,12 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes( source_regions_.scalar_flux_new(se) *= normalization_factor; } +// Normalize angular flux in the same way +#pragma omp parallel for + for (int64_t sea = 0; sea < n_source_angular_elements(); sea++) { + source_regions_.angular_flux_new(sea) *= normalization_factor; + } + // Accumulate cell-wise ray length tallies collected this iteration, then // update the simulation-averaged cell-wise volume estimates #pragma omp parallel for @@ -1852,4 +1867,29 @@ int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const return mesh_bin; } +// If tallying angular flux, this function is used to determine which angular +// bin the current ray contributes to. Rays are assigned to bins based on the +// unit-sphere Voronoi diagram generated from the "quadrature" angle set. +int FlatSourceDomain::lookup_angular_bin(Direction u) const +{ + const double x = u.x; + const double y = u.y; + const double z = u.z; + + double max_dot = -INFTY; + std::size_t best_bin = 0; + for (std::size_t j = 0; j < nangles_; ++j) { + const std::size_t k = 3 * j; + const double dot = x * angular_bin_angles_[k] + + y * angular_bin_angles_[k + 1] + + z * angular_bin_angles_[k + 2]; + if (dot > max_dot) { + max_dot = dot; + best_bin = j; + } + } + + return static_cast(best_bin); +} + } // namespace openmc diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index b4701ed1fa9..b74f6e9adf8 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -113,7 +113,7 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) } } -void LinearSourceDomain::normalize_scalar_flux_and_volumes( +void LinearSourceDomain::normalize_flux_and_volumes( double total_active_distance_per_iteration) { double normalization_factor = 1.0 / total_active_distance_per_iteration; diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index f88768984f2..a6f8592daab 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -961,7 +961,7 @@ SourceSite RandomRay::sample_s2() int RandomRay::angular_bin() { if (angular_bin_ == C_NONE) { - angular_bin_ = domain_->get_angular_bin(u()); + angular_bin_ = domain_->lookup_angular_bin(u()); } return angular_bin_; } diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 0a1ed0381d9..11bd8eec8c0 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -445,7 +445,7 @@ void RandomRaySimulation::simulate() domain_->finalize_discovered_source_regions(); // Normalize scalar flux and update volumes - domain_->normalize_scalar_flux_and_volumes( + domain_->normalize_flux_and_volumes( settings::n_particles * RandomRay::distance_active_); // Add source to scalar flux, compute number of FSR hits diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index badb9107733..0db8d7f1568 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -24,6 +24,7 @@ #include "openmc/tallies/filter_material.h" #include "openmc/tallies/filter_materialfrom.h" #include "openmc/tallies/filter_mesh.h" +#include "openmc/tallies/filter_meshangular.h" #include "openmc/tallies/filter_meshborn.h" #include "openmc/tallies/filter_meshmaterial.h" #include "openmc/tallies/filter_meshsurface.h" @@ -134,6 +135,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "mesh") { return Filter::create(id); + } else if (type == "meshangular") { + return Filter::create(id); } else if (type == "meshborn") { return Filter::create(id); } else if (type == "meshmaterial") { diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index a0698992d01..f713367bcd0 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -259,7 +259,8 @@ extern "C" int openmc_mesh_filter_get_rotation( // Check the filter type const auto& filter = model::tally_filters[index]; - if (filter->type() != FilterType::MESH) { + if (filter->type() != FilterType::MESH && + filter->type() != FilterType::MESH_ANGULAR) { set_errmsg("Tried to get a rotation from a non-mesh filter."); return OPENMC_E_INVALID_TYPE; } @@ -281,7 +282,8 @@ extern "C" int openmc_mesh_filter_set_rotation( const auto& filter = model::tally_filters[index]; // Check the filter type - if (filter->type() != FilterType::MESH) { + if (filter->type() != FilterType::MESH && + filter->type() != FilterType::MESH_ANGULAR) { set_errmsg("Tried to set a rotation from a non-mesh filter."); return OPENMC_E_INVALID_TYPE; } diff --git a/src/tallies/filter_meshangular.cpp b/src/tallies/filter_meshangular.cpp new file mode 100644 index 00000000000..c0409a9eaf6 --- /dev/null +++ b/src/tallies/filter_meshangular.cpp @@ -0,0 +1,49 @@ +#include "openmc/tallies/filter_meshsurface.h" + +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/mesh.h" + +namespace openmc { + +void MeshAngularFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + Direction u = p.u(); + if (!rotation_.empty()) { + u = u.rotate(rotation_); + } + auto bin = model::meshes[mesh_]->get_bin(r); + if (bin >= 0) { + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } +} + +void MeshAngularFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", model::meshes[mesh_]->id_); + if (rotated_) { + write_dataset(filter_group, "rotation", rotation_); + } +} + +//============================================================================== +// C-API functions +//============================================================================== + +extern "C" int openmc_meshangular_filter_get_mesh( + int32_t index, int32_t* index_mesh) +{ + return openmc_mesh_filter_get_mesh(index, index_mesh); +} + +extern "C" int openmc_meshangular_filter_set_mesh( + int32_t index, int32_t index_mesh) +{ + return openmc_mesh_filter_set_mesh(index, index_mesh); +} + +} // namespace openmc From 566415faf463bd872136e8470cb05de0cd45e8a3 Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Fri, 21 Aug 2026 18:05:38 -0400 Subject: [PATCH 4/7] Add pointset generators --- openmc/model/funcs.py | 335 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 335 insertions(+) diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index e076b080a9d..320bfbfd432 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -1,5 +1,7 @@ from collections.abc import Iterable from math import sqrt +from itertools import product +import numpy as np from operator import attrgetter from warnings import warn @@ -290,3 +292,336 @@ def pin(surfaces, items, subdivisions=None, divide_vols=True, regions = subdivide(surfaces) cells = [Cell(fill=f, region=r) for r, f in zip(regions, items)] return Universe(cells=cells, **kwargs) + +# Pointset generators for discrete-ordinate angular meshes + +def generate_levelsymmetric_sn(N: int, mu1_sq: float = None): + """ + Generates the direction vectors of an order-N level-symmetric + quadrature set. Does not provide the corresponding weights. + + Parameters + ---------- + N : int + Even quadrature order (e.g. 2, 4, 6, ..., 20). + mu1_sq : float + Square of the first direction cosine mu_1, used to generate other + levels. Must lie in the open interval (0, 1/3). For the special case + N == 2, mu1_sq must be equal to 1/3. + + Returns + ------- + np.ndarray + Array of shape (N*(N+2), 3) giving the [x, y, z] unit-vector + directions of every point in the quadrature set. + """ + if N < 2 or N % 2 != 0: + raise ValueError(f"N must be a positive even integer, got {N}") + M = N // 2 + + # lookup mu1_sq if not provided + if mu1_sq is not None: + if N == 2: + if not np.isclose(mu1_sq, 1 / 3): + raise ValueError("For N=2, mu1_sq must equal 1/3 (got {mu1_sq}).") + mu = np.array([np.sqrt(1 / 3)]) + elif N > 20: + raise ValueError( + f"Level-symmetric quadrature generation only supported for 2<=N<=20; got {N}") + else: + if not (0 < mu1_sq < 1 / 3): + raise ValueError( + f"mu1_sq={mu1_sq} is out of the valid range (0, 1/3) for N={N}." + ) + + Delta = (1 - 3 * mu1_sq) / (M - 1) + mu_sq = np.array([mu1_sq + i * Delta for i in range(M)]) + mu = np.sqrt(mu_sq) + + else: + match N: + case 2: + pass + case 4: + mu1_sq = 0.1225148226554413 + case 6: + mu1_sq = 0.0710944373419735 + case 8: + mu1_sq = 0.0476190476190470 + case 10: + mu1_sq = 0.0358425646593916 + case 12: + mu1_sq = 0.0279600712640057 + case 14: + mu1_sq = 0.0230997020840970 + case 16: + mu1_sq = 0.0193090131285642 + case 18: + mu1_sq = 0.0167300008552435 + case 20: + mu1_sq = 0.0145451663522475 + case _: + raise ValueError( + f"Level-symmetric quadrature generation only supported for 2<=N<=20; got {N}") + + if M == 1: + mu = np.array([np.sqrt(1 / 3)]) + else: + Delta = (1 - 3 * mu1_sq) / (M - 1) + mu_sq = np.array([mu1_sq + i * Delta for i in range(M)]) + mu = np.sqrt(mu_sq) + + # generate the angles for 1 octant of unit sphere + octant_points = [] + for l in range(0, M): + for m in range(0, M): + for n in range(0, M): + if l + m + n == (M + 1): + octant_points.append((mu[l], mu[m], mu[n])) + octant_points = np.array(octant_points) + + # reflect into other octants + signs = list(product([1, -1], repeat=3)) + all_points = np.array( + [pt * np.array(s) for pt in octant_points for s in signs] + ) + + return all_points + +def generate_tcl_sn(N: int): + """ + Generate the direction vectors of an order-N triangular + Chebyshev-Legendre (TCL) quadrature set. Supports even N >= 4. + + Parameters + ---------- + N : int + Even quadrature order (e.g. 4, 6, 8, ...). + + Returns + ------- + np.ndarray + Array of shape (N*(N+2), 3) giving the [x, y, z] unit-vector + directions of every point in the quadrature set. + """ + if N < 4 or N % 2 != 0: + raise ValueError(f"N must be an even integer >= 4, got {N}") + + M = N // 2 + + # get polar levels + nodes, _ = np.polynomial.legendre.leggauss(N) + mu = np.sort(nodes[nodes > 0]) + + # generate the angles for 1 octant using Chebyshev quadrature for + # azimuthal angles + octant_points = [] + for i in range(M): + count = M - i # rings get smaller towards pole + m = mu[i] + sin_theta = np.sqrt(1 - m**2) + for k in range(1, count + 1): + phi = (2 * k - 1) * (np.pi / 2) / (2 * count) + x = sin_theta * np.cos(phi) + y = sin_theta * np.sin(phi) + z = m + octant_points.append((x, y, z)) + octant_points = np.array(octant_points) + + # reflect into other octants + signs = list(product([1, -1], repeat=3)) + all_points = np.array( + [pt * np.array(s) for pt in octant_points for s in signs] + ) + + assert np.allclose(np.linalg.norm(all_points, axis=1), 1.0), ( + "Not all generated points lie on the unit sphere" + ) + + return all_points + +def _subdivide_icosahedron_faces(vertices, faces, nu): + """ + Given a list of the coordinates of the vertices of a unit icosahedron, + and a list linking sets of these vertices to individual faces of the unit + icosahedron, this function will subdivides each edge of the icosahedron + into nu equal segments, adding vertices on the edges and faces to produce + triangular subfaces of equal size. + + Parameters + ---------- + vertices : numpy array of shape (n_verts, 3) + Coordinates of each vertex of the unit icosahedron + faces : numpy array of shape (n_faces, 3) + List of vertex indices corresponding to each face of the base + icosahedron + nu : int + Subdivision frequency, integer > 1 + + Returns + ------- + subvertices : numpy array of shape (n_verts + n_faces*(nu+1)*(nu-1)/2, 3) + List of vertices on subdivided icosahedron + subfaces : numpy array of shape (n_faces*nu**2, 3) + List of vertex indices corresponding to each face of the subdivided + icosahedron + """ + edges = np.vstack([faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [0, 2]]]) + edges = np.unique(np.sort(edges, axis=1), axis=0) + + n_faces = faces.shape[0] + n_vertices = vertices.shape[0] + n_edges = edges.shape[0] + n_int_verts = (nu - 1) * (nu - 2) // 2 + + n_subverts = n_vertices + n_edges * (nu - 1) + n_faces * n_int_verts + subvertices = np.empty((n_subverts, 3)) + subvertices[:n_vertices] = vertices + + # populate edge vertices: + # position of the k-th vertex along edge AB is given by + # (1 − w_k)·a + w_k·b, w_k = (k+1)/nu + w = np.arange(1, nu) / nu + vA = vertices[edges[:, 0]] + vB = vertices[edges[:, 1]] + edge_verts = (1 - w)[:, None, None] * vA[None] + w[:, None, None] * vB[None] + subvertices[n_vertices : n_vertices + n_edges * (nu - 1)] = ( + edge_verts.transpose(1, 0, 2).reshape(-1, 3) + ) + + f_A, f_B, f_C = faces[:, 0], faces[:, 1], faces[:, 2] + + edge_dict = {(int(a), int(b)): i for i, (a, b) in enumerate(edges)} + edge_dict.update({(int(b), int(a)): ~i for i, (a, b) in enumerate(edges)}) + + def directed_edge_indices(u_arr, v_arr): + # Return (n_faces, nu-1) global vertex indices along the u->v edge + ei = np.array([edge_dict[(int(u), int(v))] for u, v in zip(u_arr, v_arr)]) + base = n_vertices + np.where(ei >= 0, ei, ~ei) * (nu - 1) + idx = base[:, None] + np.arange(nu - 1) + idx[ei < 0] = idx[ei < 0, ::-1] + return idx + + AB = directed_edge_indices(f_A, f_B) + AC = directed_edge_indices(f_A, f_C) + BC = directed_edge_indices(f_B, f_C) + + # global indices of subvertices: + # (0,0) = corner "A," (nu, 0) = corner B, (nu, nu) = corner C + # and along the edges above + local_idx = np.empty((n_faces, nu + 1, nu + 1), dtype=int) + local_idx[:, 0, 0 ] = f_A + local_idx[:, nu, 0 ] = f_B + local_idx[:, nu, nu ] = f_C + local_idx[:, 1:nu, 0 ] = AB + r_e = np.arange(1, nu) + local_idx[:, r_e, r_e ] = AC + local_idx[:, nu, 1:nu] = BC + + # row, column indices of interior points + r_int = np.array([r for r in range(2, nu) for _ in range(r - 1)]) # (n_int,) + c_int = np.array([c for r in range(2, nu) for c in range(1, r)]) + T_base = n_vertices + n_edges * (nu - 1) + + if n_int_verts > 0: + T_start = T_base + np.arange(n_faces) * n_int_verts + local_idx[:, r_int, c_int] = ( + T_start[:, None] + np.arange(n_int_verts)[None, :] + ) + + # populate subfaces with vertex coordination info: + # local vertex (r, c) with 0 ≤ c ≤ r ≤ nu has barycentric weights + # A'=(nu-r)/nu, B'=(r-c)/nu, C'=c/nu + # interior vertices have 0 < col < row < nu + tri_list = [] + for i in range(nu): + for j in range(i): + tri_list.append([(i,j), (i+1,j), (i+1,j+1)]) + tri_list.append([(i,j), (i+1,j+1), (i,j+1)]) + tri_list.append([(i,i), (i+1,i), (i+1,i+1)]) + tri_arr = np.array(tri_list) + tr, tc = tri_arr[:, :, 0], tri_arr[:, :, 1] + subfaces = local_idx[:, tr, tc].reshape(n_faces * nu ** 2, 3) + + if n_int_verts > 0: + alpha = (nu - r_int) / nu + beta = (r_int - c_int) / nu + gamma = c_int / nu + int_verts = ( alpha[None, :, None] * vertices[f_A][:, None, :] + + beta [None, :, None] * vertices[f_B][:, None, :] + + gamma[None, :, None] * vertices[f_C][:, None, :]) + subvertices[T_base:] = int_verts.reshape(-1, 3) + + return subvertices, subfaces + +def generate_icosphere_sn(nu: int = 1, point_type: str = "centroids"): + """ + Generates direction vectors from the vertices or centroids of a + unit spherical icosahedron with principal faces subdivided at the nu-th + frequency. + + That is, beginning from a "parent" icosahedron inscribed in the unit + sphere, the edges are first divided into nu equal segments, and then the + endpoints of these segments are connected to subdivide each "parent" face + into nu^2 triangular subfaces. Lastly, the vectors describing the + locations of either the vertices or centroids of these faces are + projected back onto the surface of the unit sphere. + + Parameters + ---------- + nu : int + Subdivision frequency + point_type: str + Keyword specifying whether to return the vertices ("vertices", + "vertex", or "vert") or centroids ("centroids", "centroid", or + "cent") of the subtriangles on the icosphere + + Returns + ------- + pointset : numpy array of shape (12 + 10 * (nu+1) * (nu-1), 3) if + specifying "vertices" or of shape (20 * nu**2, 3) if specifying + "centroids" + + """ + # check pointset type + match point_type: + case "vertices" | "vertex" | "vert" : + return_type = "vert" + case "centroids" | "centroid" | "cent": + return_type = "cent" + case _: + return ValueError(f"Unknown pointset type specified (got {point_type})") + + # vertices of base icosahedron + phi = (1 + np.sqrt(5)) / 2 + vertices = np.array([ + [0, 1, phi], [0, -1, phi], [1, phi, 0], + [-1, phi, 0], [phi, 0, 1], [-phi, 0, 1], + [0, -1, -phi], [0, 1, -phi], [-1, -phi, 0], + [1, -phi, 0], [-phi, 0, -1], [phi, 0, -1]]) + vertices /= np.sqrt(1 + phi ** 2) + + # coordination of vertices to faces: + # each entry in the array is a row vector, corresponding to one + # individual face of the icosahedron, giving the indices into the + # "vertices" array of its own particular vertices + faces = np.array([ + [0, 5, 1], [0, 3, 5], [0, 2, 3], [0, 4, 2], [0, 1, 4], + [1, 5, 8], [5, 3, 10], [3, 2, 7], [2, 4, 11], [4, 1, 9], + [7, 11, 6], [11, 9, 6], [9, 8, 6], [8, 10, 6], [10, 7, 6], + [2, 11, 7], [4, 9, 11], [1, 8, 9], [5, 10, 8], [3, 7, 10]]) + + # subdividing + if nu > 1: + vertices, faces = _subdivide_icosahedron_faces(vertices, faces, nu) + # project back to unit length + vertices = vertices / np.sqrt(np.sum(vertices ** 2, axis=1, keepdims=True)) + + if return_type == "vert": + pointset = vertices + else: + # return centroids + pointset = vertices[faces].mean(axis=1) + pointset = pointset / np.sqrt(np.sum(pointset ** 2, axis=1, keepdims=True)) + + return pointset \ No newline at end of file From f82335804388132fb823bf0365d909796a89fdb7 Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Tue, 25 Aug 2026 10:52:04 -0400 Subject: [PATCH 5/7] remove mesh data field --- include/openmc/mesh.h | 3 +-- openmc/mesh.py | 38 +++----------------------------------- openmc/model/funcs.py | 6 +++--- src/mesh.cpp | 6 ------ 4 files changed, 7 insertions(+), 46 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 9c10356b1d6..79aca5eb822 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -891,7 +891,7 @@ class UnitSpherePointset : public AngularMesh { void material_volumes(int nx, int ny, int nz, int max_materials, int32_t* materials, double* volumes, double* bboxes) const override { - fatal_error("material_volumes() is not supported for UnitSpherePointset") + fatal_error("material_volumes() is not supported for UnitSpherePointset"); } std::string get_mesh_type() const override { return mesh_type; } @@ -900,7 +900,6 @@ class UnitSpherePointset : public AngularMesh { void to_hdf5_inner(hid_t group) const override; vector points_; - vector data_; }; #ifdef OPENMC_DAGMC_ENABLED diff --git a/openmc/mesh.py b/openmc/mesh.py index 9d73a6a4a87..46b17630727 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -3401,8 +3401,6 @@ class UnitSpherePointset(AngularMesh): Unit-vector endpoints constituting the mesh. Must be either a numpy array or a nested list and have shape (N,3), where each row provides the [x, y, z] components of one such vector. - data : iterable of float - Scalar data (e.g. flux) assigned to each point in the mesh mesh_id : int Unique identifier for the mesh name : str @@ -3415,20 +3413,16 @@ class UnitSpherePointset(AngularMesh): name : str Name of the mesh points : numpy array of float - Array storing the - data : numpy array of float - Scalar data assigned to each point in the mesh + Array storing the direction vectors constituting the mesh. """ def __init__( self, points, - data: float | None = None, mesh_id: int | None = None, name: str = '', ): super().__init__(mesh_id, name) self.points = points - self.data = data @property def points(self): @@ -3442,20 +3436,6 @@ def points(self, pts): raise ValueError( "Unit vector array for UnitSpherePointset must have shape (N,3).") self._points = pts - - @property - def data(self): - return self._data - - @data.setter - def data(self, d): - if d is not None: - cv.check_type("AngularMesh data", d, Iterable, Real) - d = np.asarray(d).flatten() - if len(d) != len(self.points): - raise ValueError( - "Data for UnitSpherePointset does not match number of mesh elements.") - self._data = d @property def n_elements(self): @@ -3491,9 +3471,7 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): n = points.size // 3 points = points.reshape(n, 3) - data = group['data'][()] if 'data' in group else None - - return cls(points, data=data, mesh_id=mesh_id, name=name) + return cls(points, mesh_id=mesh_id, name=name) def to_xml_element(self): """Return XML representation of the mesh @@ -3513,10 +3491,6 @@ def to_xml_element(self): subelement = ET.SubElement(element, "points") subelement.text = ' '.join(map(str, pts)) - if self.data is not None: - subelement = ET.SubElement(element, "data") - subelement.text = ' '.join(map(str, self.data)) - return element @classmethod @@ -3540,13 +3514,7 @@ def from_xml_element(cls, elem: ET.Element): n = points.size // 3 points = points.reshape(n, 3) - data_elem = elem.find("data") - if data_elem is not None: - data = get_elem_list(elem, "data", float) - else: - data = None - - return cls(points, data, mesh_id=mesh_id) + return cls(points, mesh_id=mesh_id) def _read_meshes(elem): diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index 320bfbfd432..784ce5e4f77 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -295,7 +295,7 @@ def pin(surfaces, items, subdivisions=None, divide_vols=True, # Pointset generators for discrete-ordinate angular meshes -def generate_levelsymmetric_sn(N: int, mu1_sq: float = None): +def levelsymmetric_sn(N: int, mu1_sq: float = None): """ Generates the direction vectors of an order-N level-symmetric quadrature set. Does not provide the corresponding weights. @@ -388,7 +388,7 @@ def generate_levelsymmetric_sn(N: int, mu1_sq: float = None): return all_points -def generate_tcl_sn(N: int): +def tcl_sn(N: int): """ Generate the direction vectors of an order-N triangular Chebyshev-Legendre (TCL) quadrature set. Supports even N >= 4. @@ -554,7 +554,7 @@ def directed_edge_indices(u_arr, v_arr): return subvertices, subfaces -def generate_icosphere_sn(nu: int = 1, point_type: str = "centroids"): +def icosphere_sn(nu: int = 1, point_type: str = "centroids"): """ Generates direction vectors from the vertices or centroids of a unit spherical icosahedron with principal faces subdivided at the nu-th diff --git a/src/mesh.cpp b/src/mesh.cpp index e288b455967..b9d0d76d8d2 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2474,9 +2474,6 @@ UnitSpherePointset::UnitSpherePointset(hid_t group) : AngularMesh(group) points_.reserve(n); for (int i = 0; i < n; ++i) points_.push_back({flat[3 * i], flat[3 * i + 1], flat[3 * i + 2]}); - - if (object_exists(group, "data")) - read_dataset(group, "data", data_); } void UnitSpherePointset::to_hdf5_inner(hid_t mesh_group) const @@ -2489,9 +2486,6 @@ void UnitSpherePointset::to_hdf5_inner(hid_t mesh_group) const flat[3 * i + 2] = points_[i].z; } write_dataset(mesh_group, "points", flat); - - if (!data_.empty()) - write_dataset(mesh_group, "data", data_); } int UnitSpherePointset::get_bin(Direction u) const From 2cdebea935c12a81b6c289fb810ec4adeecf422e Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Wed, 26 Aug 2026 11:21:11 -0400 Subject: [PATCH 6/7] Compact data format for angular_flux_new, angular quadrature storage --- .../openmc/random_ray/flat_source_domain.h | 26 +- include/openmc/random_ray/source_region.h | 36 ++- src/random_ray/flat_source_domain.cpp | 232 +++++++++++++++--- src/random_ray/source_region.cpp | 31 ++- 4 files changed, 260 insertions(+), 65 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 216fe0ec436..9906bf98675 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -7,6 +7,7 @@ #include "openmc/random_ray/parallel_map.h" #include "openmc/random_ray/source_region.h" #include "openmc/source.h" +#include "openmc/mesh.h" #include #include @@ -36,6 +37,7 @@ class FlatSourceDomain { int64_t add_source_to_scalar_flux(); virtual void batch_reset(); void convert_source_regions_to_tallies(int64_t start_sr_id); + void initialize_angular_quadrature(); void reset_tally_volumes(); void random_ray_tally(); virtual void accumulate_iteration_flux(); @@ -68,16 +70,14 @@ class FlatSourceDomain { { return source_regions_.n_source_regions() * negroups_; } - int64_t n_source_angular_elements() const - { - return source_regions_.n_source_angular_elements(); - } 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; int lookup_mesh_idx(int64_t sr) const; - int lookup_angular_bin(Direction u) const; - Direction angular_mesh_direction(int a) { return angular_bin_angles_[a]; } + int get_angular_bin(Direction u) const; + Direction angular_quadrature_direction(int a) const; + bool tally_angular_flux_applies( + int cell_idx, int material, int mesh_idx) const; //---------------------------------------------------------------------------- // Static Data members @@ -180,12 +180,14 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Private data members int negroups_; // Number of energy groups in simulation - int nangles_; // Number of bins for any angular flux tallies - - vector angular_bin_angles_; // Directions corresponding to each - // bin midpoint in the angular flux - // tallying scheme, flattened to 1D: - // [x0, y0, z0, x1, y1, z1,...] + int nangles_ {1}; // Number of bins for any angular flux tallies + const UnitSpherePointset* angular_mesh_ {nullptr}; + + vector tally_is_angular_; + std::unordered_set angular_target_cells_; + std::unordered_set angular_target_materials_; + std::unordered_set angular_target_meshes_; + bool tally_angular_flux_everywhere_ {false}; double simulation_volume_; // Total physical volume of the simulation domain, as diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 8698f0a1427..d6726f3fa86 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -166,6 +166,7 @@ class SourceRegionHandle { double* volume_naive_; int* position_recorded_; int* external_source_present_; + int* needs_angular_flux_; Position* position_; Position* centroid_; Position* centroid_iteration_; @@ -244,6 +245,12 @@ class SourceRegionHandle { return *external_source_present_; } + int& needs_angular_flux() { return *needs_angular_flux_; } + const int needs_angular_flux() const + { + return *needs_angular_flux_; + } + Position& position() { return *position_; } const Position position() const { return *position_; } @@ -332,7 +339,8 @@ class SourceRegion { public: //---------------------------------------------------------------------------- // Constructors - SourceRegion(int negroups, bool is_linear, int nangles = 1); + SourceRegion(int negroups, bool is_linear, int nangles = 1, + bool needs_angular_flux = false); SourceRegion() = default; //---------------------------------------------------------------------------- @@ -355,6 +363,8 @@ class SourceRegion { int position_recorded_ {0}; //!< Has the position been recorded yet? int external_source_present_ { 0}; //!< Is an external source present in this region? + int needs_angular_flux_ { + 0}; //!< Is angular flux tallying active in this region? int is_small_ {0}; //!< Is it "small", receiving < 1.5 hits per iteration? int n_hits_ {0}; //!< Number of total hits (ray crossings) // Mesh that subdivides this source region @@ -470,6 +480,15 @@ class SourceRegionContainer { return external_source_present_[sr]; } + int& needs_angular_flux(int64_t sr) + { + return needs_angular_flux_[sr]; + } + const int needs_angular_flux(int64_t sr) const + { + return needs_angular_flux_[sr]; + } + Position& position(int64_t sr) { return position_[sr]; } const Position position(int64_t sr) const { return position_[sr]; } @@ -594,11 +613,11 @@ class SourceRegionContainer { double& angular_flux_new(int64_t sr, int g, int a) { - return angular_flux_new_[index_angle(sr, g, a)]; + return angular_flux_new_[angular_flux_offset_[sr] + g * nangles_ + a]; } const double angular_flux_new(int64_t sr, int g, int a) const { - return angular_flux_new_[index_angle(sr, g, a)]; + return angular_flux_new_[angular_flux_offset_[sr] + g * nangles_ + a]; } double& angular_flux_new(int64_t sea) { return angular_flux_new_[sea]; } const double angular_flux_new(int64_t sea) const @@ -660,10 +679,7 @@ class SourceRegionContainer { void flux_swap(); int64_t n_source_regions() const { return n_source_regions_; } int64_t n_source_elements() const { return n_source_regions_ * negroups_; } - int64_t n_source_angular_elements() const - { - return n_source_regions_ * negroups_ * nangles_; - } + int64_t n_source_angular_elements() const { return angular_flux_new_.size(); } int& negroups() { return negroups_; } const int negroups() const { return negroups_; } int& nangles() { return nangles_; } @@ -697,6 +713,7 @@ class SourceRegionContainer { vector volume_naive_; vector position_recorded_; vector external_source_present_; + vector needs_angular_flux_; vector position_; vector centroid_; vector centroid_iteration_; @@ -713,6 +730,7 @@ class SourceRegionContainer { vector scalar_flux_new_; vector scalar_flux_final_; vector angular_flux_new_; + vector angular_flux_offset_; vector source_; vector external_source_; @@ -733,10 +751,6 @@ class SourceRegionContainer { // Helper function for indexing inline int64_t index(int64_t sr, int g) const { return sr * negroups_ + g; } - inline int64_t index_angle(int64_t sr, int g, int a) const - { - return (sr * negroups_ + g) * nangles_ + a; - } }; } // namespace openmc diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 94376c9fefb..8ced8417220 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -5,6 +5,7 @@ #include "openmc/eigenvalue.h" #include "openmc/geometry.h" #include "openmc/material.h" +#include "openmc/mesh.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/output.h" @@ -12,6 +13,10 @@ #include "openmc/random_ray/random_ray.h" #include "openmc/simulation.h" #include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_cell.h" +#include "openmc/tallies/filter_material.h" +#include "openmc/tallies/filter_mesh.h" +#include "openmc/tallies/filter_mesh_angular.h" #include "openmc/tallies/tally.h" #include "openmc/tallies/tally_scoring.h" #include "openmc/timer.h" @@ -54,9 +59,9 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) } } - // Count the number of angular bins if performing angular flux tallying - // and store the angle set. - nangles = + // Count the number of angular bins for each SourceRegion if performing + // angular flux tallying and store the angle set. + initialize_angular_quadrature(); // Initialize source regions. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; @@ -98,7 +103,8 @@ void FlatSourceDomain::batch_reset() } #pragma omp parallel for - for (int64_t sea = 0; sea < n_source_angular_elements(); sea++) { + for (int64_t sea = 0; sea < source_regions_.n_source_angular_elements(); + sea++) { source_regions_.angular_flux_new(sea) = 0.0; } } @@ -190,7 +196,8 @@ void FlatSourceDomain::normalize_flux_and_volumes( // Normalize angular flux in the same way #pragma omp parallel for - for (int64_t sea = 0; sea < n_source_angular_elements(); sea++) { + for (int64_t sea = 0; sea < source_regions_.n_source_angular_elements(); + sea++) { source_regions_.angular_flux_new(sea) *= normalization_factor; } @@ -243,6 +250,7 @@ void FlatSourceDomain::set_flux_to_source(int64_t sr, int g) // Combine transport flux contributions and flat source contributions from the // previous iteration to generate this iteration's estimate of scalar flux. +// Not performed for angular flux as sources are isotropic int64_t FlatSourceDomain::add_source_to_scalar_flux() { int64_t n_hits = 0; @@ -404,6 +412,115 @@ void FlatSourceDomain::compute_k_eff() k_eff_ = k_eff_new; } +// Scan all tallies for a MeshAngularFilter and determine where angular flux +// tallying will occur and on how many angles +void FlatSourceDomain::initialize_angular_quadrature() +{ + angular_mesh_ = nullptr; + tally_is_angular_.assign(model::tallies.size(), false); + angular_target_cells_.clear(); + angular_target_materials_.clear(); + angular_target_meshes_.clear(); + tally_angular_flux_everywhere_ = false; + + for (int i_tally = 0; i_tally < model::tallies.size(); i_tally++) { + Tally& tally {*model::tallies[i_tally]}; + + bool has_angular_filter = false; + bool has_spatial_filter = false; + + vector tally_cells; + vector tally_materials; + vector tally_meshes; + + for (auto i_filt : tally.filters()) { + Filter* filt = model::tally_filters[i_filt].get(); + + if (auto* angular_filter = dynamic_cast(filt)) { + has_angular_filter = true; + + const auto* this_mesh = dynamic_cast( + model::meshes[angular_filter->mesh()].get()); + if (!this_mesh) { + fatal_error("MeshAngularFilter must reference a UnitSpherePointset " + "mesh for tallying in Random Ray."); + } + // if we've already found an angular tally that uses a different mesh, + // throw error + if (angular_mesh_ && angular_mesh_ != this_mesh) { + fatal_error( + "Multiple distinct angular quadratures found across " + "MeshAngularFilter tallies. All angularly-resolved tallies in " + "a Random Ray simulation must reference the same " + "UnitSpherePointset mesh."); + } + angular_mesh_ = this_mesh; + continue; + } + + if (auto* cell_filter = dynamic_cast(filt)) { + has_spatial_filter = true; + for (int32_t c : cell_filter->cells()) { + tally_cells.push_back(c); + } + } else if (auto* mat_filter = dynamic_cast(filt)) { + has_spatial_filter = true; + for (int32_t m : mat_filter->materials()) { + tally_materials.push_back(m); + } + } else if (auto* mesh_filter = dynamic_cast(filt)) { + has_spatial_filter = true; + tally_meshes.push_back(mesh_filter->mesh()); + } + } + + if (!has_angular_filter) + continue; + + tally_is_angular_[i_tally] = true; + + if (!has_spatial_filter) { + tally_angular_flux_everywhere_ = true; + continue; + } + + for (int32_t c : tally_cells) { + angular_target_cells_.insert(c); + } + for (int32_t m : tally_materials) { + angular_target_materials_.insert(m); + } + for (int32_t m : tally_meshes) { + angular_target_meshes_.insert(m); + } + } + + nangles_ = angular_mesh_ ? angular_mesh_->points_.size() : 1; +} + +// Determines whether a source region with the given cell/material/mesh +// should have angular flux tallied. +bool FlatSourceDomain::tally_angular_flux_applies( + int cell_idx, int material, int mesh_idx) const +{ + if (tally_angular_flux_everywhere_) { + return true; + } + if (angular_target_cells_.count(cell_idx)) { + return true; + } + if (material != MATERIAL_VOID && + angular_target_materials_.count(material)) { + return true; + } + if (mesh_idx != C_NONE && angular_target_meshes_.count(mesh_idx)) { + // will flag every source region this mesh applies to, whether or + // not the region actually falls in a bin + return true; + } + return false; +} + // This function is responsible for generating a mapping between random // ray flat source regions (cell instances) and tally bins. The mapping // takes the form of a "TallyTask" object, which accounts for one single @@ -449,6 +566,9 @@ void FlatSourceDomain::convert_source_regions_to_tallies(int64_t start_sr_id) // Tracks if we've generated a mapping yet for all source regions. bool all_source_regions_mapped = true; + // Get flags for tallies that need a loop layer to cover angular bins + const vector& tally_is_angular = tally_is_angular_; + // Attempt to generate mapping for all source regions #pragma omp parallel for for (int64_t sr = start_sr_id; sr < n_source_regions(); sr++) { @@ -490,6 +610,11 @@ void FlatSourceDomain::convert_source_regions_to_tallies(int64_t start_sr_id) // to what happens when scanning for applicable tallies during // MC transport. for (int i_tally = 0; i_tally < model::tallies.size(); i_tally++) { + // skip angle-dependent tallies for now + if (tally_is_angular[i_tally]) { + continue; + } + Tally& tally {*model::tallies[i_tally]}; // Initialize an iterator over valid filter bin combinations. @@ -522,6 +647,37 @@ void FlatSourceDomain::convert_source_regions_to_tallies(int64_t start_sr_id) // Reset all the filter matches for the next tally event. for (auto& match : p.filter_matches()) match.bins_present_ = false; + + // Now loop over angle-dependent tallies + if (source_regions_.needs_angular_flux(sr)) { + for (int a = 0; a < source_regions_.nangles(); a++) { + p.u() = angular_quadrature_direction(a); + + for (int i_tally = 0; i_tally < model::tallies.size(); i_tally++) { + if (!tally_is_angular[i_tally]) + continue; + Tally& tally {*model::tallies[i_tally]}; + + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true, &p.filter_matches()); + if (filter_iter == end) + continue; + + for (; filter_iter != end; ++filter_iter) { + auto filter_index = filter_iter.index_; + + for (int score = 0; score < tally.scores_.size(); score++) { + auto score_bin = tally.scores_[score]; + TallyTask task(i_tally, filter_index, score, score_bin, a); + source_regions_.tally_task(sr, g).push_back(task); + source_regions_.volume_task(sr).insert(task); + } + } + } + for (auto& match : p.filter_matches()) + match.bins_present_ = false; + } + } } } openmc::simulation::time_tallies.stop(); @@ -658,17 +814,24 @@ void FlatSourceDomain::random_ray_tally() // Determine numerical score value for (auto& task : source_regions_.tally_task(sr, g)) { + // If the current task belongs to an angle-dependent tally, replace + // scalar flux with angular flux at the task's angle. + double task_flux = (task.angle_bin == C_NONE) + ? flux + : source_regions_.angular_flux_new(sr, g, task.angle_bin) * + source_normalization_factor; + double score = 0.0; switch (task.score_type) { case SCORE_FLUX: - score = flux * volume; + score = task_flux * volume; break; case SCORE_TOTAL: if (material != MATERIAL_VOID) { score = - flux * volume * + task_flux * volume * sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] * density_mult; } @@ -677,7 +840,7 @@ void FlatSourceDomain::random_ray_tally() case SCORE_FISSION: if (material != MATERIAL_VOID) { score = - flux * volume * + task_flux * volume * sigma_f_[(material * ntemperature_ + temp) * negroups_ + g] * density_mult; } @@ -686,7 +849,7 @@ void FlatSourceDomain::random_ray_tally() case SCORE_NU_FISSION: if (material != MATERIAL_VOID) { score = - flux * volume * + task_flux * volume * nu_sigma_f_[(material * ntemperature_ + temp) * negroups_ + g] * density_mult; } @@ -698,7 +861,7 @@ void FlatSourceDomain::random_ray_tally() case SCORE_KAPPA_FISSION: score = - flux * volume * + task_flux * volume * kappa_fission_[(material * ntemperature_ + temp) * negroups_ + g] * density_mult; break; @@ -1642,13 +1805,6 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // Additionally, we need to determine the source region's material, initialize // the starting scalar flux guess, and apply any known external sources. - // Call the basic constructor for the source region and store in the parallel - // map. - bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; - SourceRegion* sr_ptr = - discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); - SourceRegionHandle handle {*sr_ptr}; - // Determine the material int gs_i_cell = gs.lowest_coord().cell(); Cell& cell = *model::cells[gs_i_cell]; @@ -1671,6 +1827,18 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } } + // Determine whether to size the source region's angular flux array for + // tallying + bool needs_angular_flux = + tally_angular_flux_applies(gs_i_cell, material, mesh_idx); + + // Call the basic constructor for the source region and store in the parallel + // map. + bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; + SourceRegion* sr_ptr = discovered_source_regions_.emplace( + sr_key, {negroups_, is_linear, nangles_, needs_angular_flux}); + SourceRegionHandle handle {*sr_ptr}; + handle.material() = material; handle.temperature_idx() = temp; @@ -1869,27 +2037,19 @@ int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const // If tallying angular flux, this function is used to determine which angular // bin the current ray contributes to. Rays are assigned to bins based on the -// unit-sphere Voronoi diagram generated from the "quadrature" angle set. -int FlatSourceDomain::lookup_angular_bin(Direction u) const +// unit-sphere Voronoi diagram generated from the "quadrature" angle set, +// stored in the referenced angular mesh. +int FlatSourceDomain::get_angular_bin(Direction u) const { - const double x = u.x; - const double y = u.y; - const double z = u.z; - - double max_dot = -INFTY; - std::size_t best_bin = 0; - for (std::size_t j = 0; j < nangles_; ++j) { - const std::size_t k = 3 * j; - const double dot = x * angular_bin_angles_[k] + - y * angular_bin_angles_[k + 1] + - z * angular_bin_angles_[k + 2]; - if (dot > max_dot) { - max_dot = dot; - best_bin = j; - } - } + return angular_mesh_->get_bin(u); +} - return static_cast(best_bin); +// Returns the representative direction of angular quadrature bin a. Used in +// convert_source_regions_to_tallies() to link angular filter bins to specific +// quadrature angles. +Direction FlatSourceDomain::angular_quadrature_direction(int a) const +{ + return angular_mesh_->points_[a]; } } // namespace openmc diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 5a305cfe27e..103a9b99a32 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -22,6 +22,7 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) volume_naive_(&sr.volume_naive_), position_recorded_(&sr.position_recorded_), external_source_present_(&sr.external_source_present_), + needs_angular_flux_(&sr.needs_angular_flux_), position_(&sr.position_), centroid_(&sr.centroid_), centroid_iteration_(&sr.centroid_iteration_), centroid_t_(&sr.centroid_t_), mom_matrix_(&sr.mom_matrix_), mom_matrix_t_(&sr.mom_matrix_t_), @@ -41,7 +42,8 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) //============================================================================== // SourceRegion implementation //============================================================================== -SourceRegion::SourceRegion(int negroups, bool is_linear, int nangles) +SourceRegion::SourceRegion(int negroups, bool is_linear, int nangles, + bool needs_angular_flux) { if (settings::run_mode == RunMode::EIGENVALUE) { // If in eigenvalue mode, set starting flux to guess of 1 @@ -56,7 +58,11 @@ SourceRegion::SourceRegion(int negroups, bool is_linear, int nangles) scalar_flux_new_.assign(negroups, 0.0); source_.assign(negroups, 0.0); scalar_flux_final_.assign(negroups, 0.0); - angular_flux_new_.assign(negroups * nangles, 0.0f); + // Only allocate angular flux storage for regions that need it + needs_angular_flux_= needs_angular_flux; + if (needs_angular_flux_) { + angular_flux_new_.assign(negroups * nangles, 0.0f); + } tally_task_.resize(negroups); if (is_linear) { @@ -89,6 +95,7 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) volume_naive_.push_back(sr.volume_naive_); position_recorded_.push_back(sr.position_recorded_); external_source_present_.push_back(sr.external_source_present_); + needs_angular_flux_.push_back(sr.needs_angular_flux_); position_.push_back(sr.position_); volume_task_.push_back(sr.volume_task_); mesh_.push_back(sr.mesh_); @@ -126,8 +133,13 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) } // Angle- and energy-dependent flux - for (int ga = 0; ga < negroups_ * nangles_; ++ga) { - angular_flux_new_.push_back(sr.angular_flux_new_[ga]); + if (!sr.angular_flux_new_.empty()) { + angular_flux_offset_.push_back(angular_flux_new_.size()); + for (int ga = 0; ga < negroups_ * nangles_; ++ga) { + angular_flux_new_.push_back(sr.angular_flux_new_[ga]); + } + } else { + angular_flux_offset_.push_back(C_NONE); } } @@ -149,6 +161,7 @@ void SourceRegionContainer::assign( volume_naive_.clear(); position_recorded_.clear(); external_source_present_.clear(); + needs_angular_flux_.clear(); position_.clear(); mesh_.clear(); parent_sr_.clear(); @@ -165,6 +178,7 @@ void SourceRegionContainer::assign( scalar_flux_new_.clear(); scalar_flux_final_.clear(); angular_flux_new_.clear(); + angular_flux_offset_.clear(); source_.clear(); external_source_.clear(); @@ -211,6 +225,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) handle.volume_naive_ = &volume_naive(sr); handle.position_recorded_ = &position_recorded(sr); handle.external_source_present_ = &external_source_present(sr); + handle.needs_angular_flux_ = &needs_angular_flux(sr); handle.position_ = &position(sr); handle.volume_task_ = &volume_task(sr); handle.mesh_ = &mesh(sr); @@ -224,7 +239,11 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) handle.external_source_ = nullptr; } handle.scalar_flux_final_ = &scalar_flux_final(sr, 0); - handle.angular_flux_new_ = &angular_flux_new(sr, 0, 0); + if (angular_flux_offset_[sr] != C_NONE) { + handle.angular_flux_new_ = &angular_flux_new(sr, 0, 0); + } else { + handle.angular_flux_new_ = nullptr; + } handle.tally_task_ = &tally_task(sr, 0); if (handle.is_linear_) { @@ -267,7 +286,7 @@ void SourceRegionContainer::adjoint_reset() std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0); } std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0); - std::fill(angular_flux_new_.begin(), angular_flux_new_.end(), 0.0f); + std::fill(angular_flux_new_.begin(), angular_flux_new_.end(), 0.0); std::fill(source_.begin(), source_.end(), 0.0f); std::fill(external_source_.begin(), external_source_.end(), 0.0f); std::fill(source_gradients_.begin(), source_gradients_.end(), From afce8f4a501962a2a9783a708a5852eb5344649c Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Wed, 26 Aug 2026 14:00:57 -0400 Subject: [PATCH 7/7] formatting and MeshAngularFilter in Python API --- include/openmc/mesh.h | 4 +- .../openmc/random_ray/flat_source_domain.h | 6 +- include/openmc/random_ray/source_region.h | 18 +-- include/openmc/tallies/filter_meshangular.h | 3 +- openmc/filter.py | 152 +++++++++++++++++- src/random_ray/flat_source_domain.cpp | 32 ++-- src/random_ray/source_region.cpp | 18 +-- 7 files changed, 185 insertions(+), 48 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 79aca5eb822..3a4a0cd8ac2 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -826,7 +826,7 @@ class UnstructuredMesh : public Mesh { virtual void initialize() = 0; }; -// Abstract class for meshes over directions (points on the unit sphere), +// Abstract class for meshes over directions (points on the unit sphere), // as opposed to meshes over volumes. class AngularMesh : public Mesh { public: @@ -873,7 +873,7 @@ class UnitSpherePointset : public AngularMesh { UnitSpherePointset(hid_t group); //! TODO: add sampling from within a spherical Voronoi cell - Position sample_element(int32_t bin, uint64_t* seed) const override + Position sample_element(int32_t bin, uint64_t* seed) const override { fatal_error( "Sampling over a UnitSpherePointset angular mesh is not supported"); diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 9906bf98675..b3b0ab6277c 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -2,12 +2,12 @@ #define OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H #include "openmc/constants.h" +#include "openmc/mesh.h" #include "openmc/openmp_interface.h" #include "openmc/position.h" #include "openmc/random_ray/parallel_map.h" #include "openmc/random_ray/source_region.h" #include "openmc/source.h" -#include "openmc/mesh.h" #include #include @@ -179,8 +179,8 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Private data members - int negroups_; // Number of energy groups in simulation - int nangles_ {1}; // Number of bins for any angular flux tallies + int negroups_; // Number of energy groups in simulation + int nangles_ {1}; // Number of bins for any angular flux tallies const UnitSpherePointset* angular_mesh_ {nullptr}; vector tally_is_angular_; diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index d6726f3fa86..e273a0e26db 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -54,7 +54,7 @@ struct TallyTask { int64_t filter_idx; int score_idx; int score_type; - // Angular quadrature bin index. Defaults to C_NONE for angle-independent, + // Angular quadrature bin index. Defaults to C_NONE for angle-independent, // scoring via the source region's scalar flux. int angle_bin {C_NONE}; TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type, @@ -246,10 +246,7 @@ class SourceRegionHandle { } int& needs_angular_flux() { return *needs_angular_flux_; } - const int needs_angular_flux() const - { - return *needs_angular_flux_; - } + const int needs_angular_flux() const { return *needs_angular_flux_; } Position& position() { return *position_; } const Position position() const { return *position_; } @@ -339,7 +336,7 @@ class SourceRegion { public: //---------------------------------------------------------------------------- // Constructors - SourceRegion(int negroups, bool is_linear, int nangles = 1, + SourceRegion(int negroups, bool is_linear, int nangles = 1, bool needs_angular_flux = false); SourceRegion() = default; @@ -362,9 +359,9 @@ class SourceRegion { double volume_naive_ {0.0}; //!< Volume as integrated from this iteration only int position_recorded_ {0}; //!< Has the position been recorded yet? int external_source_present_ { - 0}; //!< Is an external source present in this region? + 0}; //!< Is an external source present in this region? int needs_angular_flux_ { - 0}; //!< Is angular flux tallying active in this region? + 0}; //!< Is angular flux tallying active in this region? int is_small_ {0}; //!< Is it "small", receiving < 1.5 hits per iteration? int n_hits_ {0}; //!< Number of total hits (ray crossings) // Mesh that subdivides this source region @@ -480,10 +477,7 @@ class SourceRegionContainer { return external_source_present_[sr]; } - int& needs_angular_flux(int64_t sr) - { - return needs_angular_flux_[sr]; - } + int& needs_angular_flux(int64_t sr) { return needs_angular_flux_[sr]; } const int needs_angular_flux(int64_t sr) const { return needs_angular_flux_[sr]; diff --git a/include/openmc/tallies/filter_meshangular.h b/include/openmc/tallies/filter_meshangular.h index 8535f5bf699..066619485fe 100644 --- a/include/openmc/tallies/filter_meshangular.h +++ b/include/openmc/tallies/filter_meshangular.h @@ -6,7 +6,7 @@ namespace openmc { //============================================================================== -//! Indexes the direction of particle events to a mesh. +//! Indexes the direction of particle events to a mesh. //============================================================================== class MeshAngularFilter : public MeshFilter { @@ -22,7 +22,6 @@ class MeshAngularFilter : public MeshFilter { void to_statepoint(hid_t filter_group) const override; - //---------------------------------------------------------------------------- // Accessors diff --git a/openmc/filter.py b/openmc/filter.py index 00a2f31014c..607d387a698 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1349,13 +1349,157 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs): # Initialize a Pandas DataFrame from the mesh dictionary return pd.concat([df, pd.DataFrame(filter_dict)]) -class MeshAngularFilter(Filter): +class MeshAngularFilter(MeshFilter): """Bins tally events based on incident particle's direction, using an angular mesh. - + + Parameters + ---------- + mesh : openmc.MeshBase + The mesh object that events will be tallied onto + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + mesh : openmc.MeshBase + The mesh object that events will be tallied onto + translation : Iterable of float + This array specifies a vector that is used to translate (shift) + the mesh for this filter + id : int + Unique identifier for the filter + bins : list of tuple + A list of mesh indices / surfaces for each filter bin, e.g. [(1, 1, + 'x-min out'), (1, 1, 'x-min in'), ...]. Surface names use the mesh's + axis labels (e.g. r/phi/z for a cylindrical mesh). + num_bins : Integral + The number of filter bins + """ - def __init__(self): - pass + def __init__(self, mesh, filter_id=None): + self.mesh = mesh + self.id = filter_id + self._rotation = None + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tRotation', self.rotation) + return string + + @mesh.setter + def mesh(self, mesh): + cv.check_type('filter mesh', mesh, openmc.UnitSpherePointset) + self._mesh = mesh + self.bins = list(mesh.indices) + + def translation(self): + raise AttributeError( + "MeshAngularFilter instances do not permit translation.") + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'][()].decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'][()].decode() + " instead") + + if 'meshes' not in kwargs: + raise ValueError(cls.__name__ + " requires a 'meshes' keyword " + "argument.") + + mesh_id = group['bins'][()] + mesh_obj = kwargs['meshes'][mesh_id] + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + + out = cls(mesh_obj, filter_id=filter_id) + + rotation = group.get('rotation') + if rotation: + out.rotation = rotation[()] + + return out + + def get_pandas_dataframe(self, data_size, stride, **kwargs): + """Builds a Pandas DataFrame for MeshAngularFilter's bins. + + This method constructs a Pandas DataFrame object for the filter with + columns annotated by filter bin information. This is a helper method for + :meth:`Tally.get_pandas_dataframe`. + + Parameters + ---------- + data_size : int + The total number of bins in the tally corresponding to this filter + stride : int + Stride in memory for the filter + + Returns + ------- + pandas.DataFrame + A Pandas DataFrame with columns describing the mesh cell indices + corresponding to each filter bin. Column names are element index for + AngularMesh). The number of rows in the DataFrame is the same + as the total number of bins in the corresponding tally, with the + filter bin appropriately tiled to map to the corresponding tally + bins. + + See also + -------- + Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() + + """ + # Initialize dictionary to build Pandas Multi-index column + filter_dict = {} + + # Append mesh ID as outermost index of multi-index + mesh_key = f'mesh {self.mesh.id}' + + # Determine index base (0-based for angular meshes) + idx_start = 0 + + # Generate a multi-index sub-column for each axis + for label, dim_size in zip(self.mesh.axis_labels, self.mesh.dimension): + filter_dict[mesh_key, label] = _repeat_and_tile( + np.arange(idx_start, idx_start + dim_size), stride, data_size) + stride *= dim_size + + return pd.DataFrame(filter_dict) + + def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + element : lxml.etree._Element + XML element containing filter data + + """ + element = ET.Element('filter') + element.set('id', str(self.id)) + element.set('type', self.short_name.lower()) + subelement = ET.SubElement(element, 'bins') + subelement.text = str(self.mesh.id) + if self.rotation is not None: + element.set('rotation', ' '.join(map(str, self.rotation.ravel()))) + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element, **kwargs) -> MeshFilter: + mesh_id = int(get_text(elem, 'bins')) + mesh_obj = kwargs['meshes'][mesh_id] + filter_id = int(get_text(elem, "id")) + out = cls(mesh_obj, filter_id=filter_id) + + rotation = get_elem_list(elem, 'rotation', float) or [] + if rotation: + if len(rotation) == 3: + out.rotation = rotation + elif len(rotation) == 9: + out.rotation = np.array(rotation).reshape(3, 3) + return out class CollisionFilter(Filter): diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 8ced8417220..237dadad9ba 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -59,7 +59,7 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) } } - // Count the number of angular bins for each SourceRegion if performing + // Count the number of angular bins for each SourceRegion if performing // angular flux tallying and store the angle set. initialize_angular_quadrature(); @@ -412,7 +412,7 @@ void FlatSourceDomain::compute_k_eff() k_eff_ = k_eff_new; } -// Scan all tallies for a MeshAngularFilter and determine where angular flux +// Scan all tallies for a MeshAngularFilter and determine where angular flux // tallying will occur and on how many angles void FlatSourceDomain::initialize_angular_quadrature() { @@ -445,7 +445,7 @@ void FlatSourceDomain::initialize_angular_quadrature() fatal_error("MeshAngularFilter must reference a UnitSpherePointset " "mesh for tallying in Random Ray."); } - // if we've already found an angular tally that uses a different mesh, + // if we've already found an angular tally that uses a different mesh, // throw error if (angular_mesh_ && angular_mesh_ != this_mesh) { fatal_error( @@ -498,7 +498,7 @@ void FlatSourceDomain::initialize_angular_quadrature() nangles_ = angular_mesh_ ? angular_mesh_->points_.size() : 1; } -// Determines whether a source region with the given cell/material/mesh +// Determines whether a source region with the given cell/material/mesh // should have angular flux tallied. bool FlatSourceDomain::tally_angular_flux_applies( int cell_idx, int material, int mesh_idx) const @@ -509,12 +509,11 @@ bool FlatSourceDomain::tally_angular_flux_applies( if (angular_target_cells_.count(cell_idx)) { return true; } - if (material != MATERIAL_VOID && - angular_target_materials_.count(material)) { + if (material != MATERIAL_VOID && angular_target_materials_.count(material)) { return true; } if (mesh_idx != C_NONE && angular_target_meshes_.count(mesh_idx)) { - // will flag every source region this mesh applies to, whether or + // will flag every source region this mesh applies to, whether or // not the region actually falls in a bin return true; } @@ -647,7 +646,7 @@ void FlatSourceDomain::convert_source_regions_to_tallies(int64_t start_sr_id) // Reset all the filter matches for the next tally event. for (auto& match : p.filter_matches()) match.bins_present_ = false; - + // Now loop over angle-dependent tallies if (source_regions_.needs_angular_flux(sr)) { for (int a = 0; a < source_regions_.nangles(); a++) { @@ -814,12 +813,13 @@ void FlatSourceDomain::random_ray_tally() // Determine numerical score value for (auto& task : source_regions_.tally_task(sr, g)) { - // If the current task belongs to an angle-dependent tally, replace + // If the current task belongs to an angle-dependent tally, replace // scalar flux with angular flux at the task's angle. - double task_flux = (task.angle_bin == C_NONE) - ? flux - : source_regions_.angular_flux_new(sr, g, task.angle_bin) * - source_normalization_factor; + double task_flux = + (task.angle_bin == C_NONE) + ? flux + : source_regions_.angular_flux_new(sr, g, task.angle_bin) * + source_normalization_factor; double score = 0.0; switch (task.score_type) { @@ -2037,15 +2037,15 @@ int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const // If tallying angular flux, this function is used to determine which angular // bin the current ray contributes to. Rays are assigned to bins based on the -// unit-sphere Voronoi diagram generated from the "quadrature" angle set, +// unit-sphere Voronoi diagram generated from the "quadrature" angle set, // stored in the referenced angular mesh. int FlatSourceDomain::get_angular_bin(Direction u) const { return angular_mesh_->get_bin(u); } -// Returns the representative direction of angular quadrature bin a. Used in -// convert_source_regions_to_tallies() to link angular filter bins to specific +// Returns the representative direction of angular quadrature bin a. Used in +// convert_source_regions_to_tallies() to link angular filter bins to specific // quadrature angles. Direction FlatSourceDomain::angular_quadrature_direction(int a) const { diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 103a9b99a32..1a975e966e0 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -22,12 +22,12 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) volume_naive_(&sr.volume_naive_), position_recorded_(&sr.position_recorded_), external_source_present_(&sr.external_source_present_), - needs_angular_flux_(&sr.needs_angular_flux_), - position_(&sr.position_), centroid_(&sr.centroid_), - centroid_iteration_(&sr.centroid_iteration_), centroid_t_(&sr.centroid_t_), - mom_matrix_(&sr.mom_matrix_), mom_matrix_t_(&sr.mom_matrix_t_), - volume_task_(&sr.volume_task_), mesh_(&sr.mesh_), - parent_sr_(&sr.parent_sr_), scalar_flux_old_(sr.scalar_flux_old_.data()), + needs_angular_flux_(&sr.needs_angular_flux_), position_(&sr.position_), + centroid_(&sr.centroid_), centroid_iteration_(&sr.centroid_iteration_), + centroid_t_(&sr.centroid_t_), mom_matrix_(&sr.mom_matrix_), + mom_matrix_t_(&sr.mom_matrix_t_), volume_task_(&sr.volume_task_), + mesh_(&sr.mesh_), parent_sr_(&sr.parent_sr_), + scalar_flux_old_(sr.scalar_flux_old_.data()), scalar_flux_new_(sr.scalar_flux_new_.data()), source_(sr.source_.data()), external_source_(sr.external_source_.data()), scalar_flux_final_(sr.scalar_flux_final_.data()), @@ -42,8 +42,8 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) //============================================================================== // SourceRegion implementation //============================================================================== -SourceRegion::SourceRegion(int negroups, bool is_linear, int nangles, - bool needs_angular_flux) +SourceRegion::SourceRegion( + int negroups, bool is_linear, int nangles, bool needs_angular_flux) { if (settings::run_mode == RunMode::EIGENVALUE) { // If in eigenvalue mode, set starting flux to guess of 1 @@ -59,7 +59,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear, int nangles, source_.assign(negroups, 0.0); scalar_flux_final_.assign(negroups, 0.0); // Only allocate angular flux storage for regions that need it - needs_angular_flux_= needs_angular_flux; + needs_angular_flux_ = needs_angular_flux; if (needs_angular_flux_) { angular_flux_new_.assign(negroups * nangles, 0.0f); }