Homec4science
Diffusion wireles-gpu (master) Browse

wireles-gpu/src/filter_gpu.cuf
8183685fbe4dmaster

Last Change
Raw File

filter_gpu.cuf
View Options

module filter_gpu
use iso_c_binding
use precision
use dimen
contains
subroutine filter_ddxy_gpu(f,dfdx,dfdy,t,plan)
!---------------------------------------------------------------------------
! declaration
!---------------------------------------------------------------------------
use cufft_m
use cudafor
implicit none
! Input: velocity field
real(fp), dimension(nx,ny,nz2),intent(inout),device :: f,dfdx, dfdy
complex(fp), dimension(:,:,:),allocatable,device :: f_hat,dfdx_hat,dfdy_hat
complex(fp), dimension(nx,ny,nz2) :: f_hat_host
! fft plan
type(c_ptr) :: plan(2)
integer :: istat,k,t
! Cuda variable
type(dim3) :: grid, tBlock
tBlock = dim3(16,16,2)
grid = dim3((nx/2)/16,ny/16,nzb/2)
save f_hat,dfdx_hat,dfdy_hat
!---------------------------------------------------------------------------
! main code
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! init
if (t.eq.1) then
allocate(f_hat(nx/2+1,ny,nz2))
allocate(dfdx_hat(nx/2+1,ny,nz2))
allocate(dfdy_hat(nx/2+1,ny,nz2))
endif
!---------------------------------------------------------------------------
! Forward fft
call cufftExec(plan(1),f(:,:,2:nzb+1),f_hat(:,:,2:nzb+1))
!---------------------------------------------------------------------------
! Apply sharp spectral filter
call filter_kernel<<<grid, tBlock>>>(f_hat,dfdx_hat,dfdy_hat)
!---------------------------------------------------------------------------
! Backward fft
call cufftExec(plan(2),f_hat(:,:,2:nzb+1),f(:,:,2:nzb+1))
call cufftExec(plan(2),dfdx_hat(:,:,2:nzb+1),dfdx(:,:,2:nzb+1))
call cufftExec(plan(2),dfdy_hat(:,:,2:nzb+1),dfdy(:,:,2:nzb+1))
end subroutine filter_ddxy_gpu
! subroutine filter_2laa (f_hat,f_hatd,f,t,plan)
! !---------------------------------------------------------------------------
! ! declaration
! !---------------------------------------------------------------------------
! use cufft_m
! use cudafor
! implicit none
! integer i,j,k,ii,jj,t,flag
! type(c_ptr),intent(in) :: plan(2)
! real(fp),dimension(nx,ny,nz2),intent(in),device :: f_hat,f_hatd,f
! complex(fp), dimension(:,:,:),allocatable,device :: f_hat_ft,f_hat_ftd
! type(dim3) :: grid, tBlock
! tBlock = dim3(16,16,2)
! grid = dim3((nx/2)/16,ny/16,nz2/2)
! save f_hat_ft,f_hat_ftd
! !---------------------------------------------------------------------------
! ! main code
! !---------------------------------------------------------------------------
! !---------------------------------------------------------------------------
! ! init
! if (t == 1) then
! allocate(f_hat_ft (nx/2,ny,nz2))
! allocate(f_hat_ftd(nx/2,ny,nz2))
! end if
! !---------------------------------------------------------------------------
! ! compute
! call cufftExec(plan(1),f(:,:,2:nzb+1),f_hat_ft(:,:,2:nzb+1))
! call filter_kernel_test<<<grid,tBlock>>>(f_hat_ft,f_hat_ftd)
! call cufftExec(plan(2),f_hat_ft (:,:,2:nzb+1),f_hat (:,:,2:nzb+1))
! call cufftExec(plan(2),f_hat_ftd(:,:,2:nzb+1),f_hatd(:,:,2:nzb+1))
! end subroutine filter_2laa
subroutine filter_2dsl(f_filtered,f,t,plan_2d)
!---------------------------------------------------------------------------
! declaration
!---------------------------------------------------------------------------
use cudafor
use cufft_m
implicit none
integer :: i,j,k,ii,jj,t,flag,iret,istat
type(c_ptr) :: plan_2d(2)
real(fp),dimension(nx,ny),intent(in),device :: f
real(fp),dimension(nx,ny),intent(out),device :: f_filtered
complex(fp), dimension(:,:),allocatable,device :: f_hat
type(dim3) :: grid, tBlock
grid= dim3((nx/2)/16,ny/16,1)
tBlock= dim3(16,16,1)
save f_hat
!---------------------------------------------------------------------------
! main code
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! init
if (t.eq.1) then
allocate(f_hat(nx/2+1,ny))
end if
!---------------------------------------------------------------------------
! compute
! forward fft --> f
call cufftExec(plan_2d(1),f,f_hat)
! Apply sharp spectral cutoff. Note: square cutoff
call filter_kernel_2d<<<grid, tBlock>>>(f_hat)
! backward fft --> f
call cufftExec(plan_2d(2),f_hat,f_filtered)
end subroutine filter_2dsl
attributes(global) subroutine filter_kernel(f_hat,dfdx_hat,dfdy_hat)
implicit none
complex(fp), dimension(nx/2+1,ny,nz2),intent(inout) :: f_hat
complex(fp), dimension(nx/2+1,ny,nz2),intent(inout) :: dfdx_hat,dfdy_hat
integer :: i,j,k,ii,jj
i = (blockIdx%x - 1) * blockDim%x + threadIdx%x
j = (blockIdx%y - 1) * blockDim%y + threadIdx%y
k = (blockIdx%z - 1) * blockDim%z + threadIdx%z
k = k + 1 ! k = 2,nzb + 1
!---------------------------------------------------------------------------
! Scale back the fft
f_hat(i,j,k) = f_hat(i,j,k)*inxny
!---------------------------------------------------------------------------
! Preparing sharp spectral cut-off
! Shift wave number
ii = i - 1
jj = j - 1
if (jj > int(ny/2)) then
jj = jj-ny
end if
jj = jj*l_r
!---------------------------------------------------------------------------
! Apply sharp spectral cut-off
call spec_box_cutoff(f_hat,fgr,i,j,k)
!---------------------------------------------------------------------------
! Apply sharp spectral cut-off
! dfdx_hat
dfdx_hat(i,j,k) = dcmplx(aimag(f_hat(i,j,k))*(-1.0), &
dble(f_hat(i,j,k)))*ii
! dfdy_hat
dfdy_hat(i,j,k) = dcmplx(aimag(f_hat(i,j,k))*(-1.0), &
dble(f_hat(i,j,k)))*jj
! nx/2 + 1 is not cover by the grid, thus needs to be treated individually
if (i==nx/2) then
dfdx_hat(nx/2+1,j,k) = dcmplx(0._fp)
dfdy_hat(nx/2+1,j,k) = dcmplx(0._fp)
endif
end subroutine filter_kernel
attributes(global) subroutine filter_kernel_2d(f_hat)
! 2d Filter at the wall
implicit none
complex(fp), dimension(nx/2+1,ny),intent(inout) :: f_hat
integer :: i,j,ii,jj
i = (blockIdx%x - 1) * blockDim%x + threadIdx%x;
j = (blockIdx%y - 1) * blockDim%y + threadIdx%y;
! Scale back the fft
f_hat(i,j) = f_hat(i,j)*inxny
! ! Shift wave number
ii = i - 1
jj = j - 1
if(jj.gt.int(nx/2)) jj=jj-ny
jj=jj*l_r
if (ii >= floor(nx/(2*fgr*tfr)))then
f_hat(i,j)= dcmplx(0.0_fp)
! filter y
elseif(abs(jj) >= floor(l_r*ny/(2*fgr*tfr)))then
f_hat(i,j)= dcmplx(0.0_fp)
end if
f_hat(nx/2+1,j) = dcmplx(0.0_fp)
end subroutine filter_kernel_2d
attributes(global) subroutine test_filter_kernel(f_hat_ft,filterRadius)
implicit none
complex(fp), dimension(nx/2+1,ny,nz2),intent(inout) :: f_hat_ft
real(fp),intent(in),value :: filterRadius
integer :: i,j,k,ii,jj
i = (blockIdx%x - 1) * blockDim%x + threadIdx%x
j = (blockIdx%y - 1) * blockDim%y + threadIdx%y
k = (blockIdx%z - 1) * blockDim%z + threadIdx%z
! Scale back the fft
f_hat_ft(i,j,k) = f_hat_ft(i,j,k)*inxny
call spec_box_cutoff(f_hat_ft,filterRadius,i,j,k)
end subroutine test_filter_kernel
attributes(device) subroutine spec_box_cutoff(f_hat,filter,i,j,k)
implicit none
complex(fp), dimension(nx/2+1,ny,nz2),intent(inout) :: f_hat
real(fp),intent(in) :: filter
integer,intent(in) :: i,j,k
integer :: ii,jj
!---------------------------------------------------------------------------
! Preparing sharp spectral cut-off
! Shift wave number
ii = i - 1
jj = j - 1
if(jj.gt.int(ny/2)) jj = jj-ny
jj = jj*l_r
!---------------------------------------------------------------------------
! Apply sharp spectral cut-off
! filter x
if (ii >= floor(nx/(2.0*filter))) then
f_hat(i,j,k)=dcmplx(0._fp)
! filter y
elseif(abs(jj) >= floor(l_r*ny/(2.0*filter))) then
f_hat(i,j,k)=dcmplx(0._fp)
end if
! nx/2 + 1 is not cover by the grid, thus needs to be treated individually
if (i==nx/2) then
f_hat(nx/2+1,j,k) = dcmplx(0._fp)
end if
end subroutine spec_box_cutoff
end module filter_gpu
c4science · Help