From fa3fc916816fb79c19f95749aec12fdf32e41f53 Mon Sep 17 00:00:00 2001 From: j-fletcher Date: Tue, 16 Jun 2026 16:05:43 -0400 Subject: [PATCH 1/4] 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/4] 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/4] 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/4] 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