Skip to content

Instantly share code, notes, and snippets.

@eschnett
Created Nov 12, 2020
Embed
What would you like to do?
Non-local boundary conditions
! Apply various non-local boundary conditions to a 3d grid
module slab
implicit none
integer, parameter :: nghosts = 2
contains
! Given an input array that is defined in the interior (on the dots
! and stars), fill the boundaries by rotating the array by 90
! degrees in either direction. The pivot point is the top left dot.
! The array values marked with a star end up being used to fill in
! the boundaries, the array values marked with a dot end up unused.
!
! Formula: f(-x,y,z) = f(y,x,z)
! f(x,-y,z) = f(y,x,z)
!
! Array shape:
! 1 2 2 2 2
! 3 . . . .
! 3 . * * *
! 3 . * * *
! 3 . * * *
subroutine rotating90(dst, src)
double precision, intent(inout) :: dst(:,:,:)
double precision, intent(in) :: src(:,:,:)
integer :: i, j, k
if (size(dst) /= size(src)) stop
if (size(dst,1) /= size(dst,2)) stop
i0 = nghosts
j0 = nghosts
! Fill region marked "1" by rotating array by 180 degrees:
forall (k = 1 : size(a,3))
forall (j = 1 : nghosts)
forall (i = 1 : nghosts)
dst(i0-i, j0-j, k) = src(i0+i, j0+j, k)
end forall
end forall
end forall
! Fill region marked "2" by rotating array 90 degrees counterclockwise:
forall (k = 1 : size(a,3))
forall (j = 1 : nghosts)
forall (i = 1 : size(a,1)-nghosts-1)
dst(i0+i, j0-j, k) = src(i0+j, j0+i, k)
end forall
end forall
end forall
! Fill region marked "3" by rotating array 90 degrees clockwise:
forall (k = 1 : size(a,3))
forall (j = 1 : size(a,2)-nghosts-1)
forall (i = 1 : nghosts)
dst(i0-i, j0+j, k) = src(i0+j, j0+i, k)
end forall
end forall
end forall
end subroutine rotating90
! Given an input array that is defined in the interior (on the dots
! and stars), fill the boundaries by rotating the array by 180
! degrees in either direction. The pivot point is the centre dot.
! The array values marked with a star end up being used to fill in
! the boundaries, the array values marked with a dot end up unused.
!
! Formula: f(-x,y,z) = f(x,-y,z)
!
! Array shape:
! 1 . * * *
! 1 . * * *
! 1 . * * *
! 1 . * * *
! 1 . * * *
subroutine rotating180(dst, src)
double precision, intent(inout) :: dst(:,:,:)
double precision, intent(in) :: src(:,:,:)
integer :: i, j, k
if (size(dst) /= size(src)) stop
if (mod(size(dst,2), 2) /= 1) stop
i0 = nghosts
j0 = nghosts
! Fill region marked "1" by rotating array by 180 degrees:
forall (k = 1 : size(a,3))
forall (j = 1 : size(a,2))
forall (i = 1 : nghosts)
dst(i0-i, j, k) = src(i0+i, size(a,2)+1-j, k)
end forall
end forall
end forall
end subroutine rotating180
! Given an input array that is defined in the interior (on the
! stars), fill the boundaries by cycling the array by 1/2 its length
! in the x direction and wrapping around in the y direction.
!
! Formula: f(-theta,phi,z) = f(theta,phi+pi,z)
!
! Array shape:
! 1 1 1 2 2 2
! . . . . . .
! * * * * * *
! * * * * * *
! * * * * * *
! . . . . . .
! 3 3 3 4 4 4
subroutine polar
double precision, intent(inout) :: dst(:,:,:)
double precision, intent(in) :: src(:,:,:)
integer :: i, j, k
if (size(dst) /= size(src)) stop
if (mod(size(dst,1), 2) /= 0) stop
i0 = nghosts
j0 = nghosts
! Fill regions marked "3" and "4":
forall (k = 1 : size(a,3))
forall (j = 1 : nghosts)
forall (i = 1 : size(a,1))
ii = mod(i + size(a,1) - 1, size(a,1)) + 1
jj = j + size(a,2) - 2*nghosts
dst(i, j0, k) = src(1 + mod(ii, jj, k)
end forall
end forall
end forall
! ditto for regions 3 and 4
end subroutine polar
end module slab
@WeiqunZhang

This comment has been minimized.

Copy link

@WeiqunZhang WeiqunZhang commented Nov 13, 2020

Should Line 47 forall (i = 1 : size(a,1)-nghosts-1) read forall (i = 1 : size(a,1)-nghosts)?

@eschnett

This comment has been minimized.

Copy link
Owner Author

@eschnett eschnett commented Nov 14, 2020

Yes, probably. Please treat this as pseudo-code to illustrate what the boundary conditions do. I didn't test the code thoroughly, and there are likely errors.

@WeiqunZhang

This comment has been minimized.

Copy link

@WeiqunZhang WeiqunZhang commented Nov 14, 2020

For rotating90, I assume you do want the ghost cells at the higher end of the domain being filled with rotated data too. Is that correct?

@eschnett

This comment has been minimized.

Copy link
Owner Author

@eschnett eschnett commented Nov 15, 2020

The boundary conditions should be applied at the lower i and j boundaries only, i.e. on 2 of the 4 boundaries.

The 180 degree rotating boundaries should be applied at the lower i boundary only.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment