diff --git a/openmc/surface.py b/openmc/surface.py index 239b868aa2..ca2d5d6b45 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -222,6 +222,9 @@ class Plane(Surface): boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. + periodic_surface : openmc.Surface + If a periodic boundary condition is used, the surface with which this + one is periodic with coefficients : dict Dictionary of surface coefficients id : int @@ -239,6 +242,7 @@ class Plane(Surface): self._type = 'plane' self._coeff_keys = ['A', 'B', 'C', 'D'] + self._periodic_surface = None self.a = A self.b = B self.c = C @@ -260,6 +264,10 @@ class Plane(Surface): def d(self): return self.coefficients['D'] + @property + def periodic_surface(self): + return self._periodic_surface + @a.setter def a(self, A): check_type('A coefficient', A, Real) @@ -280,6 +288,21 @@ class Plane(Surface): check_type('D coefficient', D, Real) self._coefficients['D'] = D + @periodic_surface.setter + def periodic_surface(self, periodic_surface): + check_type('periodic surface', periodic_surface, Plane) + self._periodic_surface = periodic_surface + periodic_surface._periodic_surface = self + + def create_xml_subelement(self): + element = super(Plane, self).create_xml_subelement() + + # Add periodic surface pair information + if self.boundary_type == 'periodic': + if self.periodic_surface is not None: + element.set("periodic_surface_id", str(self.periodic_surface.id)) + return element + class XPlane(Plane): """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0` @@ -306,6 +329,9 @@ class XPlane(Plane): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} Boundary condition that defines the behavior for particles hitting the surface. + periodic_surface : openmc.Surface + If a periodic boundary condition is used, the surface with which this + one is periodic with coefficients : dict Dictionary of surface coefficients id : int @@ -391,6 +417,9 @@ class YPlane(Plane): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} Boundary condition that defines the behavior for particles hitting the surface. + periodic_surface : openmc.Surface + If a periodic boundary condition is used, the surface with which this + one is periodic with coefficients : dict Dictionary of surface coefficients id : int @@ -477,6 +506,9 @@ class ZPlane(Plane): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} Boundary condition that defines the behavior for particles hitting the surface. + periodic_surface : openmc.Surface + If a periodic boundary condition is used, the surface with which this + one is periodic with coefficients : dict Dictionary of surface coefficients id : int diff --git a/src/geometry.F90 b/src/geometry.F90 index 2059ea42fc..21c8baa6b8 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -502,26 +502,26 @@ contains select type (surf) type is (SurfaceXPlane) - select type (opposite => surfaces(surf % opposite) % obj) + select type (opposite => surfaces(surf % i_periodic) % obj) type is (SurfaceXPlane) p % coord(1) % xyz(1) = opposite % x0 end select type is (SurfaceYPlane) - select type (opposite => surfaces(surf % opposite) % obj) + select type (opposite => surfaces(surf % i_periodic) % obj) type is (SurfaceYPlane) p % coord(1) % xyz(2) = opposite % y0 end select type is (SurfaceZPlane) - select type (opposite => surfaces(surf % opposite) % obj) + select type (opposite => surfaces(surf % i_periodic) % obj) type is (SurfaceZPlane) p % coord(1) % xyz(3) = opposite % z0 end select end select ! Reassign particle's surface - p % surface = sign(surfaces(surf % opposite) % obj % id, p % surface) + p % surface = sign(surf % i_periodic, p % surface) ! Figure out what cell particle is in now p % n_coord = 1 diff --git a/src/initialize.F90 b/src/initialize.F90 index 74cab87c4b..09bedb1388 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -580,20 +580,6 @@ contains class(Lattice), pointer :: lat => null() type(TallyObject), pointer :: t => null() - ! Adjust opposite surfaces for periodic boundaries - do i = 1, size(surfaces) - associate (surf => surfaces(i) % obj) - if (surf % bc == BC_PERIODIC) then - if (surface_dict % has_key(surf % opposite)) then - surf % opposite = surface_dict % get_key(surf % opposite) - else - call fatal_error("Could not find opposite surface " // & - trim(to_str(surf % opposite)) // ".") - end if - end if - end associate - end do - do i = 1, n_cells ! ======================================================================= ! ADJUST REGION SPECIFICATION FOR EACH CELL diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 242d3459a6..8a9b2299fb 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1578,6 +1578,12 @@ contains case ('periodic') s%bc = BC_PERIODIC boundary_exists = .true. + + ! Check for specification of periodic surface + if (check_for_node(node_surf, "periodic_surface_id")) then + call get_node_value(node_surf, "periodic_surface_id", & + s % i_periodic) + end if case default call fatal_error("Unknown boundary condition '" // trim(word) // & &"' specified on surface " // trim(to_str(s%id))) @@ -1597,33 +1603,45 @@ contains if (surfaces(i) % obj % bc == BC_PERIODIC) then select type (surf => surfaces(i) % obj) type is (SurfaceXPlane) - if (i == i_xmin) then - surf % opposite = i_xmax - elseif (i == i_xmax) then - surf % opposite = i_xmin + if (surf % i_periodic == NONE) then + if (i == i_xmin) then + surf % i_periodic = i_xmax + elseif (i == i_xmax) then + surf % i_periodic = i_xmin + else + call fatal_error("Periodic boundary condition applied to & + &interior surface.") + end if else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") + surf % i_periodic = surface_dict % get_key(surf % i_periodic) end if type is (SurfaceYPlane) - if (i == i_ymin) then - surf % opposite = i_ymax - elseif (i == i_ymax) then - surf % opposite = i_ymin + if (surf % i_periodic == NONE) then + if (i == i_ymin) then + surf % i_periodic = i_ymax + elseif (i == i_ymax) then + surf % i_periodic = i_ymin + else + call fatal_error("Periodic boundary condition applied to & + &interior surface.") + end if else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") + surf % i_periodic = surface_dict % get_key(surf % i_periodic) end if type is (SurfaceZPlane) - if (i == i_zmin) then - surf % opposite = i_zmax - elseif (i == i_zmax) then - surf % opposite = i_zmin + if (surf % i_periodic == NONE) then + if (i == i_zmin) then + surf % i_periodic = i_zmax + elseif (i == i_zmax) then + surf % i_periodic = i_zmin + else + call fatal_error("Periodic boundary condition applied to & + &interior surface.") + end if else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") + surf % i_periodic = surface_dict % get_key(surf % i_periodic) end if class default @@ -1633,7 +1651,7 @@ contains ! Make sure opposite surface is also periodic associate (surf => surfaces(i) % obj) - if (surfaces(surf % opposite) % obj % bc /= BC_PERIODIC) then + if (surfaces(surf % i_periodic) % obj % bc /= BC_PERIODIC) then call fatal_error("Could not find matching surface for periodic & &boundary on surface " // trim(to_str(surf % id)) // ".") end if diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 0b5d3c86be..68e5144b7c 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -1,6 +1,6 @@ module surface_header - use constants, only: ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT + use constants, only: NONE, ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT implicit none @@ -15,8 +15,8 @@ module surface_header neighbor_pos(:), & ! List of cells on positive side neighbor_neg(:) ! List of cells on negative side integer :: bc ! Boundary condition - integer :: opposite ! Opposite surface for periodic boundary - character(len=104) :: name = "" ! User-defined name + integer :: i_periodic = NONE ! Index of corresponding periodic surface + character(len=104) :: name = "" ! User-defined name contains procedure :: sense procedure :: reflect diff --git a/tests/test_periodic/geometry.xml b/tests/test_periodic/geometry.xml deleted file mode 100644 index 6ecfec1972..0000000000 --- a/tests/test_periodic/geometry.xml +++ /dev/null @@ -1,12 +0,0 @@ - - - - - - - - - - - - diff --git a/tests/test_periodic/inputs_true.dat b/tests/test_periodic/inputs_true.dat new file mode 100644 index 0000000000..d50d0b8592 --- /dev/null +++ b/tests/test_periodic/inputs_true.dat @@ -0,0 +1 @@ +af589996f2930337afe34ba9894098ff5efe3b29b6e927117220b718bf29b630ffdbc931754d465a8e8100125a8aa997dbe10aab322b43f69d59710573996a6d \ No newline at end of file diff --git a/tests/test_periodic/materials.xml b/tests/test_periodic/materials.xml deleted file mode 100644 index a7bf4faf4e..0000000000 --- a/tests/test_periodic/materials.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - - - - - - - - - - diff --git a/tests/test_periodic/results_true.dat b/tests/test_periodic/results_true.dat index f65dbafd1b..b0bdb22c7d 100644 --- a/tests/test_periodic/results_true.dat +++ b/tests/test_periodic/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.542742E+00 4.410461E-02 +1.040109E+00 6.527490E-02 diff --git a/tests/test_periodic/settings.xml b/tests/test_periodic/settings.xml deleted file mode 100644 index af09407ae5..0000000000 --- a/tests/test_periodic/settings.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - 1000 - 4 - 0 - - - - -5. -5. -5. 5. 5. 5. - - - diff --git a/tests/test_periodic/test_periodic.py b/tests/test_periodic/test_periodic.py index b584632f08..558514575a 100644 --- a/tests/test_periodic/test_periodic.py +++ b/tests/test_periodic/test_periodic.py @@ -3,9 +3,58 @@ import os import sys sys.path.insert(0, os.pardir) -from testing_harness import TestHarness +from testing_harness import PyAPITestHarness +import openmc + + +class PeriodicTest(PyAPITestHarness): + def _build_inputs(self): + # Define materials + water = openmc.Material(1) + water.add_nuclide('H-1', 2.0) + water.add_nuclide('O-16', 1.0) + water.add_s_alpha_beta('HH2O', '71t') + water.set_density('g/cc', 1.0) + + fuel = openmc.Material(2) + fuel.add_nuclide('U-235', 1.0) + fuel.set_density('g/cc', 4.5) + + materials = openmc.Materials((water, fuel)) + materials.default_xs = '71c' + materials.export_to_xml() + + # Define geometry + x_min = openmc.XPlane(1, x0=-5., boundary_type='periodic') + x_max = openmc.XPlane(2, x0=5., boundary_type='periodic') + x_max.periodic_surface = x_min + + y_min = openmc.YPlane(3, y0=-5., boundary_type='periodic') + y_max = openmc.YPlane(4, y0=5., boundary_type='periodic') + + z_min = openmc.ZPlane(5, z0=-5., boundary_type='reflective') + z_max = openmc.ZPlane(6, z0=5., boundary_type='reflective') + z_cyl = openmc.ZCylinder(7, x0=-2.5, y0=2.5, R=2.0) + + outside_cyl = openmc.Cell(1, fill=water, region=( + +x_min & -x_max & +y_min & -y_max & +z_min & -z_max & +z_cyl)) + inside_cyl = openmc.Cell(2, fill=fuel, region=+z_min & -z_max & -z_cyl) + root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl)) + + geometry = openmc.Geometry() + geometry.root_universe = root_universe + geometry.export_to_xml() + + # Define settings + settings = openmc.Settings() + settings.particles = 1000 + settings.batches = 4 + settings.inactive = 0 + settings.source = openmc.Source(space=openmc.stats.Box( + *outside_cyl.region.bounding_box)) + settings.export_to_xml() if __name__ == '__main__': - harness = TestHarness('statepoint.4.h5') + harness = PeriodicTest('statepoint.4.h5') harness.main()