{{ message }}

Instantly share code, notes, and snippets.

# eschnett/slab.f90

Created Nov 12, 2020
Non-local boundary conditions
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 ! 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 commented Nov 13, 2020

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

### 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 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 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.
to join this conversation on GitHub. Already have an account? Sign in to comment