crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #880 from paulromano/rotational-periodic

Browse source 
Add support for rotationally periodic BCs
This commit is contained in:
sam 2017-05-22 09:03:23 -04:00 • committed by GitHub
commit 57d58d7de8
2 changed files with 56 additions and 11 deletions
Download patch file Download diff file Expand all files Collapse all files

19
docs/source/usersguide/geometry.rst
View file

@ -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:

48
src/geometry.F90
View file

@ -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