Pair periodic surfaces for RectangularParallelepiped

openmc/model/surface_composite.py

@@ -19,18 +19,39 @@
 class CompositeSurface(ABC):
     """Multiple primitive surfaces combined into a composite surface"""
 
+    def _periodic_pairs(self):
+        """Names of component surfaces that are periodic with each other"""
+        components = {id(getattr(self, name)): name
+                      for name in self._surface_names}
+        pairs = []
+        for name in self._surface_names:
+            partner = getattr(getattr(self, name), 'periodic_surface', None)
+            if partner is not None and id(partner) in components:
+                pairs.append((name, components[id(partner)]))
+        return pairs
+
+    def _restore_periodic_pairs(self, pairs):
+        # Transforming a component can replace it with a new surface (or a
+        # copy of its partner), so point each pair back at its new partner.
+        for name, partner_name in pairs:
+            getattr(self, name).periodic_surface = getattr(self, partner_name)
+
     def translate(self, vector, inplace=False):
+        pairs = self._periodic_pairs()
         surf = self if inplace else copy(self)
         for name in self._surface_names:
             s = getattr(surf, name)
             setattr(surf, name, s.translate(vector, inplace))
+        surf._restore_periodic_pairs(pairs)
         return surf
 
     def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
+        pairs = self._periodic_pairs()
         surf = self if inplace else copy(self)
         for name in self._surface_names:
             s = getattr(surf, name)
             setattr(surf, name, s.rotate(rotation, pivot, order, inplace))
+        surf._restore_periodic_pairs(pairs)
         return surf
 
     @property
@@ -646,6 +667,19 @@ def __init__(self, xmin, xmax, ymin, ymax, zmin, zmax, **kwargs):
         self.ymax = openmc.YPlane(y0=ymax, **kwargs)
         self.zmin = openmc.ZPlane(z0=zmin, **kwargs)
         self.zmax = openmc.ZPlane(z0=zmax, **kwargs)
+        if kwargs.get('boundary_type') == 'periodic':
+            self.xmin.periodic_surface = self.xmax
+            self.ymin.periodic_surface = self.ymax
+            self.zmin.periodic_surface = self.zmax
+
+    @CompositeSurface.boundary_type.setter
+    def boundary_type(self, boundary_type):
+        for name in self._surface_names:
+            getattr(self, name).boundary_type = boundary_type
+        if boundary_type == 'periodic':
+            self.xmin.periodic_surface = self.xmax
+            self.ymin.periodic_surface = self.ymax
+            self.zmin.periodic_surface = self.zmax
 
     def __neg__(self):
         return -self.xmax & +self.xmin & -self.ymax & +self.ymin & -self.zmax & +self.zmin

tests/unit_tests/test_periodic_bc.py

@@ -45,3 +45,27 @@ def test_rotational_periodic_bc(angle):
 
     # Run the model
     model.run()
+
+
+def test_rectangular_parallelepiped_periodic(run_in_tmpdir):
+    """Regression test for issue #3966."""
+    mat = openmc.Material(density=1.0, density_units='g/cm3')
+    mat.add_nuclide('H1', 1.0)
+
+    inner = openmc.model.RectangularParallelepiped(
+        -3., 3., -3., 3., -3., 3., boundary_type='periodic')
+    outer = openmc.model.RectangularParallelepiped(
+        -4., 4., -4., 4., -4., 4., boundary_type='periodic')
+
+    cell_inner = openmc.Cell(fill=mat, region=-inner)
+    cell_outer = openmc.Cell(fill=mat, region=+inner & -outer)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell_inner, cell_outer])
+    model.settings.particles = 100
+    model.settings.batches = 2
+    model.settings.inactive = 1
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.IndependentSource(
+        energy=openmc.stats.delta_function(1.0e4)
+    )
+    model.run()

tests/unit_tests/test_surface_composite.py

@@ -54,6 +54,32 @@ def test_rectangular_parallelepiped():
     repr(s)
 
 
+def test_rectangular_parallelepiped_periodic():
+    s = openmc.model.RectangularParallelepiped(
+        -3., 3., -3., 3., -3., 3., boundary_type='periodic')
+    assert s.boundary_type == 'periodic'
+    assert s.xmin.periodic_surface == s.xmax
+    assert s.ymin.periodic_surface == s.ymax
+    assert s.zmin.periodic_surface == s.zmax
+
+    s.boundary_type = 'vacuum'
+    assert s.xmin.boundary_type == 'vacuum'
+
+    s.boundary_type = 'periodic'
+    assert s.xmin.periodic_surface == s.xmax
+
+
+@pytest.mark.parametrize("inplace", [False, True])
+def test_rectangular_parallelepiped_periodic_transform(inplace):
+    s = openmc.model.RectangularParallelepiped(
+        -3., 3., -3., 3., -3., 3., boundary_type='periodic')
+    for t in (s.translate((1., 2., 3.), inplace=inplace),
+              s.rotate((0., 0., 30.), inplace=inplace)):
+        for lo, hi in (('xmin', 'xmax'), ('ymin', 'ymax'), ('zmin', 'zmax')):
+            assert getattr(t, lo).periodic_surface is getattr(t, hi)
+            assert getattr(t, hi).periodic_surface is getattr(t, lo)
+
+
 @pytest.mark.parametrize(
     "axis, indices", [
         ("X", [0, 1, 2]),