From f43eb58c644ab766577ac479632c8cd0fb545630 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 18 May 2017 15:16:52 -0500 Subject: [PATCH 1/2] Support rotational periodic boundary conditions --- src/geometry.F90 | 48 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 9 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 5276edb24..f80bd7eb7 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -388,14 +388,15 @@ contains type(Particle), intent(inout) :: p integer, intent(in) :: last_cell ! last cell particle was in - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: norm ! "norm" of surface normal - real(8) :: d ! distance between point and plane - real(8) :: xyz(3) ! Saved global coordinate - integer :: i_surface ! index in surfaces - logical :: found ! particle found in universe? + real(8) :: u ! x-component of direction + real(8) :: v ! y-component of direction + real(8) :: w ! z-component of direction + real(8) :: norm ! "norm" of surface normal + real(8) :: d ! distance between point and plane + real(8) :: xyz(3) ! Saved global coordinate + integer :: i_surface ! index in surfaces + logical :: rotational ! if rotational periodic BC applied + logical :: found ! particle found in universe? class(Surface), pointer :: surf i_surface = abs(p % surface) @@ -515,17 +516,42 @@ contains p % coord(1) % xyz = xyz end if + rotational = .false. select type (surf) type is (SurfaceXPlane) select type (opposite => surfaces(surf % i_periodic) % obj) type is (SurfaceXPlane) p % coord(1) % xyz(1) = opposite % x0 + type is (SurfaceYPlane) + rotational = .true. + + ! Rotate direction + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + p % coord(1) % uvw(1) = v + p % coord(1) % uvw(2) = -u + + ! Rotate position + p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0 + p % coord(1) % xyz(2) = opposite % y0 end select type is (SurfaceYPlane) select type (opposite => surfaces(surf % i_periodic) % obj) type is (SurfaceYPlane) p % coord(1) % xyz(2) = opposite % y0 + type is (SurfaceXPlane) + rotational = .true. + + ! Rotate direction + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + p % coord(1) % uvw(1) = -v + p % coord(1) % uvw(2) = u + + ! Rotate position + p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0 + p % coord(1) % xyz(1) = opposite % x0 end select type is (SurfaceZPlane) @@ -550,7 +576,11 @@ contains end select ! Reassign particle's surface - p % surface = sign(surf % i_periodic, p % surface) + if (rotational) then + p % surface = surf % i_periodic + else + p % surface = sign(surf % i_periodic, p % surface) + end if ! Figure out what cell particle is in now p % n_coord = 1 From 266f2f85834cdd5cec78a8a94044c5da1d2cdc14 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 18 May 2017 15:24:01 -0500 Subject: [PATCH 2/2] Update user's guide section on boundary conditions --- docs/source/usersguide/geometry.rst | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index 18b9e7945..0c18e56ae 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -148,8 +148,23 @@ surface. To specify a vacuum boundary condition, simply change the Reflective and periodic boundary conditions can be set with the strings 'reflective' and 'periodic'. Vacuum and reflective boundary conditions can be -applied to any type of surface. Periodic boundary conditions can only be applied -to pairs of axis-aligned planar surfaces. +applied to any type of surface. Periodic boundary conditions can be applied to +pairs of planar surfaces. For axis-aligned planes, matching periodic surfaces +can be determined automatically. For non-axis-aligned planes, it is necessary to +specify pairs explicitly using the :attr:`Surface.periodic_surface` attribute as +in the following example:: + + p1 = openmc.Plane(A=0.3, B=5.0, D=1.0, boundary_type='periodic') + p2 = openmc.Plane(A=0.3, B=5.0, D=-1.0, boundary_type='periodic') + p1.periodic_surface = p2 + +Rotationally-periodic boundary conditions can be specified for a pair of +:class:`XPlane` and :class:`YPlane`; in that case, the +:attr:`Surface.periodic_surface` attribute must be specified manually as well. + +.. caution:: When using rotationally-periodic boundary conditions, your geometry + must be defined in the first quadrant, i.e., above the y-plane and + to the right of the x-plane. .. _usersguide_cells: