diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h index 2f48b6d4d48..2da714d8bd2 100644 --- a/include/openmc/tallies/filter_mesh.h +++ b/include/openmc/tallies/filter_mesh.h @@ -51,6 +51,12 @@ class MeshFilter : public Filter { virtual bool translated() const { return translated_; } + virtual void set_rotation(const vector& rotation); + + virtual const vector& rotation() const { return rotation_; } + + virtual bool rotated() const { return rotated_; } + protected: //---------------------------------------------------------------------------- // Data members @@ -58,6 +64,8 @@ class MeshFilter : public Filter { int32_t mesh_; //!< Index of the mesh bool translated_ {false}; //!< Whether or not the filter is translated Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation + bool rotated_ {false}; //!< Whether or not the filter is rotated + vector rotation_; //!< Filter rotation }; } // namespace openmc diff --git a/openmc/filter.py b/openmc/filter.py index b5926af5ad9..42227c9f327 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -807,6 +807,25 @@ class MeshFilter(Filter): translation : Iterable of float This array specifies a vector that is used to translate (shift) the mesh for this filter + rotation : Iterable of float + This array specifies the angles in degrees about the x, y, and z axes + that the filled universe should be rotated. The rotation applied + is an intrinsic rotation with specified Tait-Bryan angles. That is to + say, if the angles are :math:`(\phi, \theta, \psi)`, then the rotation + matrix applied is :math:`R_z(\psi) R_y(\theta) R_x(\phi)` or + + .. math:: + + \left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\phi \sin\psi + + \sin\phi \sin\theta \cos\psi & \sin\phi \sin\psi + \cos\phi + \sin\theta \cos\psi \\ \cos\theta \sin\psi & \cos\phi \cos\psi + + \sin\phi \sin\theta \sin\psi & -\sin\phi \cos\psi + \cos\phi + \sin\theta \sin\psi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi + \cos\theta \end{array} \right ] + + A rotation matrix can also be specified directly by setting this + attribute to a nested list (or 2D numpy array) that specifies each + element of the matrix. bins : list of tuple A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1), ...] @@ -819,6 +838,7 @@ def __init__(self, mesh, filter_id=None): self.mesh = mesh self.id = filter_id self._translation = None + self._rotation = None def __hash__(self): string = type(self).__name__ + '\n' @@ -830,6 +850,7 @@ def __repr__(self): string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) string += '{: <16}=\t{}\n'.format('\tID', self.id) string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation) + string += '{: <16}=\t{}\n'.format('\tRotation', self.rotation) return string @classmethod @@ -853,6 +874,10 @@ def from_hdf5(cls, group, **kwargs): if translation: out.translation = translation[()] + rotation = group.get('rotation') + if rotation: + out.rotation = rotation[()] + return out @property @@ -887,6 +912,15 @@ def translation(self, t): cv.check_length('mesh filter translation', t, 3) self._translation = np.asarray(t) + @property + def rotation(self): + return self._rotation + + @rotation.setter + def rotation(self, rotation): + cv.check_length('mesh filter rotation', rotation, 3) + self._rotation = np.asarray(rotation) + def can_merge(self, other): # Mesh filters cannot have more than one bin return False @@ -969,6 +1003,8 @@ def to_xml_element(self): element[0].text = str(self.mesh.id) if self.translation is not None: element.set('translation', ' '.join(map(str, self.translation))) + if self.rotation is not None: + element.set('rotation', ' '.join(map(str, self.rotation.ravel()))) return element @classmethod @@ -981,6 +1017,11 @@ def from_xml_element(cls, elem: ET.Element, **kwargs) -> MeshFilter: translation = elem.get('translation') if translation: out.translation = [float(x) for x in translation.split()] + rotation = elem.get('rotation') + if rotation: + values = [float(x) for x in rotation.split()] + if len(values) == 9: + out.rotation = np.array(values).reshape(3, 3) return out diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index b30f5e66282..b532de1c283 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -96,6 +96,14 @@ _dll.openmc_mesh_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)] _dll.openmc_mesh_filter_set_translation.restype = c_int _dll.openmc_mesh_filter_set_translation.errcheck = _error_handler +_dll.openmc_mesh_filter_get_rotation.argtypes = [c_int32, POINTER(c_double), + POINTER(c_size_t)] +_dll.openmc_mesh_filter_get_rotation.restype = c_int +_dll.openmc_mesh_filter_get_rotation.errcheck = _error_handler +_dll.openmc_mesh_filter_set_rotation.argtypes = [ + 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_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 @@ -392,6 +400,10 @@ class MeshFilter(Filter): Mesh used for the filter translation : Iterable of float 3-D coordinates of the translation vector + rotation : Iterable of float + The rotation matrix or angles of the filter mesh. This can either be + a fully specified 3 x 3 rotation matrix or an Iterable of length 3 + with the angles in degrees about the x, y, and z axes, respectively. """ filter_type = 'mesh' @@ -421,6 +433,34 @@ def translation(self): def translation(self, translation): _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation)) + @property + def rotation(self): + rotation_data = np.zeros(12) + rot_size = c_size_t() + + _dll.openmc_mesh_filter_get_rotation( + self._index, rotation_data.ctypes.data_as(POINTER(c_double)), + rot_size) + rot_size = rot_size.value + + if rot_size == 9: + return rotation_data[:rot_size].shape(3, 3) + elif rot_size in (0, 12): + # If size is 0, rotation_data[9:] will be zeros. This indicates no + # rotation and is the most straightforward way to always return + # an iterable of floats + return rotation_data[9:] + else: + raise ValueError( + f'Invalid size of rotation matrix: {rot_size}') + + @rotation.setter + def rotation(self, rotation_data): + flat_rotation = np.asarray(rotation_data, dtype=float).flatten() + + _dll.openmc_mesh_filter_set_rotation( + self._index, flat_rotation.ctypes.data_as(POINTER(c_double)), + c_size_t(len(flat_rotation))) class MeshBornFilter(Filter): """MeshBorn filter stored internally. diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 03f7da97847..d1ab02f4819 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -7,6 +7,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/mesh.h" +#include "openmc/position.h" #include "openmc/xml_interface.h" namespace openmc { @@ -31,6 +32,10 @@ void MeshFilter::from_xml(pugi::xml_node node) if (check_for_node(node, "translation")) { set_translation(get_node_array(node, "translation")); } + // Read the rotation transform. + if (check_for_node(node, "rotation")) { + set_rotation(get_node_array(node, "rotation")); + } } void MeshFilter::get_all_bins( @@ -46,6 +51,12 @@ void MeshFilter::get_all_bins( last_r -= translation(); r -= translation(); } + // apply rotation if present + if (!rotation_.empty()) { + last_r = last_r.rotate(rotation_); + r = r.rotate(rotation_); + u = u.rotate(rotation_); + } if (estimator != TallyEstimator::TRACKLENGTH) { auto bin = model::meshes[mesh_]->get_bin(r); @@ -66,6 +77,9 @@ void MeshFilter::to_statepoint(hid_t filter_group) const if (translated_) { write_dataset(filter_group, "translation", translation_); } + if (rotated_) { + write_dataset(filter_group, "rotation", rotation_); + } } std::string MeshFilter::text_label(int bin) const @@ -94,6 +108,40 @@ void MeshFilter::set_translation(const double translation[3]) this->set_translation({translation[0], translation[1], translation[2]}); } +void MeshFilter::set_rotation(const vector& rot) +{ + rotated_ = true; + + // Compute and store the rotation matrix. + rotation_.clear(); + rotation_.reserve(rot.size() == 9 ? 9 : 12); + if (rot.size() == 3) { + double phi = -rot[0] * PI / 180.0; + double theta = -rot[1] * PI / 180.0; + double psi = -rot[2] * PI / 180.0; + rotation_.push_back(std::cos(theta) * std::cos(psi)); + rotation_.push_back(-std::cos(phi) * std::sin(psi) + + std::sin(phi) * std::sin(theta) * std::cos(psi)); + rotation_.push_back(std::sin(phi) * std::sin(psi) + + std::cos(phi) * std::sin(theta) * std::cos(psi)); + rotation_.push_back(std::cos(theta) * std::sin(psi)); + rotation_.push_back(std::cos(phi) * std::cos(psi) + + std::sin(phi) * std::sin(theta) * std::sin(psi)); + rotation_.push_back(-std::sin(phi) * std::cos(psi) + + std::cos(phi) * std::sin(theta) * std::sin(psi)); + rotation_.push_back(-std::sin(theta)); + rotation_.push_back(std::sin(phi) * std::cos(theta)); + rotation_.push_back(std::cos(phi) * std::cos(theta)); + + // When user specifies angles, write them at end of vector + rotation_.push_back(rot[0]); + rotation_.push_back(rot[1]); + rotation_.push_back(rot[2]); + } else { + std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_)); + } +} + //============================================================================== // C-API functions //============================================================================== @@ -202,4 +250,48 @@ extern "C" int openmc_mesh_filter_set_translation( return 0; } +//! Return the rotation matrix of a mesh filter +extern "C" int openmc_mesh_filter_get_rotation( + int32_t index, double rot[], size_t* n) +{ + // Make sure this is a valid index to an allocated filter + if (int err = verify_filter(index)) + return err; + + // Check the filter type + const auto& filter = model::tally_filters[index]; + if (filter->type() != FilterType::MESH) { + set_errmsg("Tried to get a rotation from a non-mesh filter."); + return OPENMC_E_INVALID_TYPE; + } + // Get rotation from the mesh filter and set value + auto mesh_filter = dynamic_cast(filter.get()); + *n = mesh_filter->rotation().size(); + std::memcpy(rot, mesh_filter->rotation().data(), + *n * sizeof(mesh_filter->rotation()[0])); + return 0; +} + +//! Set the flattened rotation matrix of a mesh filter +extern "C" int openmc_mesh_filter_set_rotation( + int32_t index, const double rot[], size_t rot_len) +{ + // Make sure this is a valid index to an allocated filter + if (int err = verify_filter(index)) + return err; + + const auto& filter = model::tally_filters[index]; + // Check the filter type + if (filter->type() != FilterType::MESH) { + set_errmsg("Tried to set a rotation from a non-mesh filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Get a pointer to the filter and downcast + auto mesh_filter = dynamic_cast(filter.get()); + std::vector vec_rot(rot, rot + rot_len); + mesh_filter->set_rotation(vec_rot); + return 0; +} + } // namespace openmc diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 68f7550ae5a..b62e8a96556 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -818,7 +818,8 @@ void WeightWindowsGenerator::create_tally() for (const auto& f : model::tally_filters) { if (f->type() == FilterType::MESH) { const auto* mesh_filter = dynamic_cast(f.get()); - if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated()) { + if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated() && + !mesh_filter->rotated()) { ww_tally->add_filter(f.get()); found_mesh_filter = true; break; diff --git a/tests/regression_tests/filter_rotations/__init__.py b/tests/regression_tests/filter_rotations/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/filter_rotations/inputs_true.dat b/tests/regression_tests/filter_rotations/inputs_true.dat new file mode 100644 index 00000000000..5ba0017149d --- /dev/null +++ b/tests/regression_tests/filter_rotations/inputs_true.dat @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 0 + + + + 3 4 5 + -9 -9 -9 + 9 9 9 + + + 3 4 5 + -9 -9 -9 + 9 9 9 + + + 1 + + + 2 + + + 1 + total + + + 2 + total + + + diff --git a/tests/regression_tests/filter_rotations/results_true.dat b/tests/regression_tests/filter_rotations/results_true.dat new file mode 100644 index 00000000000..ba7047368ff --- /dev/null +++ b/tests/regression_tests/filter_rotations/results_true.dat @@ -0,0 +1,244 @@ +k-combined: +2.298294E-01 3.256961E-01 +tally 1: +4.880089E-02 +4.894233E-04 +8.221192E-02 +1.373712E-03 +5.205261E-02 +5.584673E-04 +1.272758E-01 +3.371623E-03 +3.599375E-01 +2.655645E-02 +1.377572E-01 +3.835778E-03 +1.646205E-01 +5.916595E-03 +3.806765E-01 +2.947576E-02 +1.470573E-01 +4.360385E-03 +5.568294E-02 +6.368003E-04 +6.775790E-02 +9.447074E-04 +5.580909E-02 +6.303584E-04 +6.674847E-02 +9.080853E-04 +9.747628E-02 +1.934325E-03 +6.824062E-02 +9.908844E-04 +1.891692E-01 +7.211874E-03 +5.976566E-01 +7.206211E-02 +1.805640E-01 +6.599680E-03 +2.114692E-01 +9.202185E-03 +6.514433E-01 +8.523087E-02 +1.905740E-01 +7.287607E-03 +7.633486E-02 +1.216488E-03 +1.117747E-01 +2.538705E-03 +7.438042E-02 +1.127109E-03 +7.356007E-02 +1.090617E-03 +1.161258E-01 +2.734941E-03 +6.745417E-02 +9.371031E-04 +2.153323E-01 +9.293477E-03 +9.414362E-01 +2.261845E-01 +2.016933E-01 +8.270477E-03 +2.370758E-01 +1.150628E-02 +9.973982E-01 +2.489516E-01 +2.152505E-01 +9.424546E-03 +7.502684E-02 +1.208304E-03 +1.244515E-01 +3.168682E-03 +7.515812E-02 +1.151397E-03 +6.355178E-02 +8.423316E-04 +1.040915E-01 +2.224064E-03 +6.902234E-02 +9.843766E-04 +1.978591E-01 +7.894161E-03 +5.162186E-01 +5.357742E-02 +1.897326E-01 +7.231935E-03 +1.827950E-01 +6.839560E-03 +5.780661E-01 +6.799910E-02 +1.881215E-01 +7.093746E-03 +6.083700E-02 +7.861240E-04 +9.978991E-02 +2.022733E-03 +6.004653E-02 +7.428556E-04 +4.951220E-02 +4.952385E-04 +7.401227E-02 +1.124154E-03 +5.151090E-02 +5.339359E-04 +1.356118E-01 +3.711749E-03 +3.400140E-01 +2.367518E-02 +1.385305E-01 +4.060592E-03 +1.438711E-01 +4.191824E-03 +3.276170E-01 +2.210924E-02 +1.279501E-01 +3.435254E-03 +4.852204E-02 +4.919417E-04 +7.229090E-02 +1.062894E-03 +4.344414E-02 +3.825354E-04 +tally 2: +4.852204E-02 +4.919417E-04 +7.229090E-02 +1.062894E-03 +4.344414E-02 +3.825354E-04 +1.438711E-01 +4.191824E-03 +3.276170E-01 +2.210924E-02 +1.279501E-01 +3.435254E-03 +1.356118E-01 +3.711749E-03 +3.400140E-01 +2.367518E-02 +1.385305E-01 +4.060592E-03 +4.951220E-02 +4.952385E-04 +7.401227E-02 +1.124154E-03 +5.151090E-02 +5.339359E-04 +6.083700E-02 +7.861240E-04 +9.978991E-02 +2.022733E-03 +6.004653E-02 +7.428556E-04 +1.827950E-01 +6.839560E-03 +5.780661E-01 +6.799910E-02 +1.881215E-01 +7.093746E-03 +1.978591E-01 +7.894161E-03 +5.162186E-01 +5.357742E-02 +1.897326E-01 +7.231935E-03 +6.355178E-02 +8.423316E-04 +1.040915E-01 +2.224064E-03 +6.902234E-02 +9.843766E-04 +7.502684E-02 +1.208304E-03 +1.244515E-01 +3.168682E-03 +7.515812E-02 +1.151397E-03 +2.370758E-01 +1.150628E-02 +9.973982E-01 +2.489516E-01 +2.152505E-01 +9.424546E-03 +2.153323E-01 +9.293477E-03 +9.414362E-01 +2.261845E-01 +2.016933E-01 +8.270477E-03 +7.356007E-02 +1.090617E-03 +1.161258E-01 +2.734941E-03 +6.745417E-02 +9.371031E-04 +7.633486E-02 +1.216488E-03 +1.117747E-01 +2.538705E-03 +7.438042E-02 +1.127109E-03 +2.114692E-01 +9.202185E-03 +6.514433E-01 +8.523087E-02 +1.905740E-01 +7.287607E-03 +1.891692E-01 +7.211874E-03 +5.976566E-01 +7.206211E-02 +1.805640E-01 +6.599680E-03 +6.674847E-02 +9.080853E-04 +9.747628E-02 +1.934325E-03 +6.824062E-02 +9.908844E-04 +5.568294E-02 +6.368003E-04 +6.775790E-02 +9.447074E-04 +5.580909E-02 +6.303584E-04 +1.646205E-01 +5.916595E-03 +3.806765E-01 +2.947576E-02 +1.470573E-01 +4.360385E-03 +1.272758E-01 +3.371623E-03 +3.599375E-01 +2.655645E-02 +1.377572E-01 +3.835778E-03 +4.880089E-02 +4.894233E-04 +8.221192E-02 +1.373712E-03 +5.205261E-02 +5.584673E-04 diff --git a/tests/regression_tests/filter_rotations/test.py b/tests/regression_tests/filter_rotations/test.py new file mode 100644 index 00000000000..e1fda44ec89 --- /dev/null +++ b/tests/regression_tests/filter_rotations/test.py @@ -0,0 +1,73 @@ +import numpy as np + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + + model = openmc.model.Model() + + fuel = openmc.Material() + fuel.set_density('g/cm3', 10.0) + fuel.add_nuclide('U235', 1.0) + zr = openmc.Material() + zr.set_density('g/cm3', 1.0) + zr.add_nuclide('Zr90', 1.0) + model.materials.extend([fuel, zr]) + + box1 = openmc.model.RectangularPrism(10.0, 10.0) + box2 = openmc.model.RectangularPrism(20.0, 20.0, boundary_type='reflective') + top = openmc.ZPlane(z0=10.0, boundary_type='vacuum') + bottom = openmc.ZPlane(z0=-10.0, boundary_type='vacuum') + cell1 = openmc.Cell(fill=fuel, region=-box1 & +bottom & -top) + cell2 = openmc.Cell(fill=zr, region=+box1 & -box2 & +bottom & -top) + model.geometry = openmc.Geometry([cell1, cell2]) + + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 1000 + + rotation = np.array((180, 0, 0)) + + llc = np.array([-9, -9, -9]) + urc = np.array([9, 9, 9]) + + mesh_dims = (3, 4, 5) + + filters = [] + + # un-rotated meshes + reg_mesh = openmc.RegularMesh() + reg_mesh.dimension = mesh_dims + reg_mesh.lower_left = llc + reg_mesh.upper_right = urc + + filters.append(openmc.MeshFilter(reg_mesh)) + + + # rotated meshes + rotated_reg_mesh = openmc.RegularMesh() + rotated_reg_mesh.dimension = mesh_dims + rotated_reg_mesh.lower_left = llc + rotated_reg_mesh.upper_right = urc + + filters.append(openmc.MeshFilter(rotated_reg_mesh)) + filters[-1].rotation = rotation + + # Create tallies + for f in filters: + tally = openmc.Tally() + tally.filters = [f] + tally.scores = ['total'] + model.tallies.append(tally) + + return model + + +def test_filter_mesh_rotations(model): + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 64c16c238ea..1394b2ba553 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -583,6 +583,13 @@ def test_regular_mesh(lib_init): assert isinstance(mesh, openmc.lib.RegularMesh) assert mesh_id == mesh.id + rotation = (180.0, 0.0, 0.0) + + mf = openmc.lib.MeshFilter(mesh) + assert mf.mesh == mesh + mf.rotation = rotation + assert np.allclose(mf.rotation, rotation) + translation = (1.0, 2.0, 3.0) mf = openmc.lib.MeshFilter(mesh)