diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3d11d39017a..3a4a0cd8ac2 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -826,6 +826,82 @@ 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_; +}; + #ifdef OPENMC_DAGMC_ENABLED class MOABMesh : public UnstructuredMesh { 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/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 09414fd4465..b3b0ab6277c 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -2,6 +2,7 @@ #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" @@ -30,12 +31,13 @@ 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(); 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(); @@ -72,6 +74,10 @@ class FlatSourceDomain { 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 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 @@ -173,7 +179,15 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Private data members - int negroups_; // Number of energy groups in simulation + 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_; + 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/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 b61d2d67aa8..e4d0bb0d443 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 1d2bbe1e8dc..e273a0e26db 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -54,9 +54,13 @@ struct TallyTask { int64_t filter_idx; int score_idx; int score_type; - TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type) + // 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, + 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 +69,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 +81,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 +147,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}; @@ -159,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_; @@ -180,6 +188,7 @@ class SourceRegionHandle { float* source_; float* external_source_; double* scalar_flux_final_; + double* angular_flux_new_; //!< only accumulated for ext. source biasing MomentArray* source_gradients_; MomentArray* flux_moments_old_; @@ -236,6 +245,9 @@ 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_; } @@ -279,6 +291,15 @@ 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 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]; } @@ -315,7 +336,8 @@ class SourceRegion { public: //---------------------------------------------------------------------------- // Constructors - SourceRegion(int negroups, bool is_linear); + SourceRegion(int negroups, bool is_linear, int nangles = 1, + bool needs_angular_flux = false); SourceRegion() = default; //---------------------------------------------------------------------------- @@ -337,7 +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? 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 @@ -375,6 +399,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 @@ -395,8 +422,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; @@ -450,6 +477,12 @@ 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]; } @@ -572,6 +605,20 @@ class SourceRegionContainer { return scalar_flux_final_[se]; } + double& angular_flux_new(int64_t sr, int g, int 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_[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 + { + return angular_flux_new_[sea]; + } + 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]; } @@ -626,8 +673,11 @@ 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 angular_flux_new_.size(); } 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); @@ -638,6 +688,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) @@ -656,6 +707,7 @@ class SourceRegionContainer { vector volume_naive_; vector position_recorded_; vector external_source_present_; + vector needs_angular_flux_; vector position_; vector centroid_; vector centroid_iteration_; @@ -671,6 +723,8 @@ class SourceRegionContainer { vector scalar_flux_old_; vector scalar_flux_new_; vector scalar_flux_final_; + vector angular_flux_new_; + vector angular_flux_offset_; vector source_; vector external_source_; 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..066619485fe --- /dev/null +++ b/include/openmc/tallies/filter_meshangular.h @@ -0,0 +1,40 @@ +#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..607d387a698 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,158 @@ 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(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, 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): """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..46b17630727 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,166 @@ 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. + 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 direction vectors constituting the mesh. + """ + def __init__( + self, + points, + mesh_id: int | None = None, + name: str = '', + ): + super().__init__(mesh_id, name) + self.points = points + + @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 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) + + return cls(points, 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)) + + 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) + + return cls(points, mesh_id=mesh_id) + def _read_meshes(elem): """Generate dictionary of meshes from a given XML node diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index e076b080a9d..784ce5e4f77 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 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 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 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 diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index d166f124795..3e3a49c32cf 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -42,6 +42,7 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const // normalize reflected u to ensure no floating point error leads to // unnormalized directions u /= u.norm(); + p.direction_changed(); // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); @@ -60,6 +61,7 @@ void WhiteBC::handle_particle(Particle& p, const Surface& surf) const // normalize outgoing u to ensure no floating point error leads to // unnormalized directions u /= u.norm(); + p.direction_changed(); // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); @@ -246,6 +248,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(); // normalize new_u to ensure no floating point error leads to unnormalized // directions diff --git a/src/mesh.cpp b/src/mesh.cpp index a0e497613c0..b9d0d76d8d2 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2429,6 +2429,80 @@ 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]}); +} + +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); +} + +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..237dadad9ba 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,13 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) } } + // 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; - source_regions_ = SourceRegionContainer(negroups_, is_linear); + source_regions_ = SourceRegionContainer(negroups_, is_linear, nangles_); // Initialize tally volumes if (volume_normalized_flux_tallies_) { @@ -92,6 +101,12 @@ 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 < source_regions_.n_source_angular_elements(); + sea++) { + source_regions_.angular_flux_new(sea) = 0.0; + } } void FlatSourceDomain::accumulate_iteration_flux() @@ -165,7 +180,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 +194,13 @@ 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 < source_regions_.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 @@ -228,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; @@ -389,6 +412,114 @@ 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 @@ -434,6 +565,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++) { @@ -475,6 +609,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. @@ -507,6 +646,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(); @@ -643,17 +813,25 @@ 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; } @@ -662,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; } @@ -671,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; } @@ -683,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; @@ -1627,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]; @@ -1656,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; @@ -1852,4 +2035,21 @@ 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, +// 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 +// quadrature angles. +Direction FlatSourceDomain::angular_quadrature_direction(int a) const +{ + return angular_mesh_->points_[a]; +} + } // 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 dde5023e44f..a6f8592daab 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_->lookup_angular_bin(u()); + } + return angular_bin_; +} + } // namespace openmc 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/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 78543c5ab53..1a975e966e0 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -10,22 +10,28 @@ 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_), - 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()), + 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()), @@ -36,7 +42,8 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) //============================================================================== // SourceRegion implementation //============================================================================== -SourceRegion::SourceRegion(int negroups, bool is_linear) +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 @@ -51,6 +58,11 @@ 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); + // 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) { @@ -83,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_); @@ -118,6 +131,16 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) // Tally tasks tally_task_.emplace_back(sr.tally_task_[g]); } + + // Angle- and energy-dependent flux + 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); + } } void SourceRegionContainer::assign( @@ -138,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(); @@ -153,6 +177,8 @@ void SourceRegionContainer::assign( scalar_flux_old_.clear(); scalar_flux_new_.clear(); scalar_flux_final_.clear(); + angular_flux_new_.clear(); + angular_flux_offset_.clear(); source_.clear(); external_source_.clear(); @@ -184,6 +210,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); @@ -198,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); @@ -211,6 +239,11 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) handle.external_source_ = nullptr; } handle.scalar_flux_final_ = &scalar_flux_final(sr, 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_) { @@ -253,6 +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.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(), 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