Homec4science
Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (3D_parallel) Browse

gyacomo/src/fft_mod.F90
9a072e2a62053D_parallel

Last Change
Raw File

fft_mod.F90
View Options

module fft
! Fast fourier transform module
! Used to supress high frequency instabilities
use prec_const
use array
implicit none
private
! Public Functions
PUBLIC :: fft_recursive, initialize_filter, suppress_high_freq
contains
! In place Cooley-Tukey FFT, tsource : rosettacode.org/wiki/Fast_Fourier_transform
recursive subroutine fft_recursive(x)
complex(kind=dp), dimension(:), intent(inout) :: x
complex(kind=dp) :: t
integer :: N
integer :: i
complex(kind=dp), dimension(:), allocatable :: even, odd
N=size(x)
if(N .le. 1) return
allocate(odd((N+1)/2))
allocate(even(N/2))
! divide
odd =x(1:N:2)
even=x(2:N:2)
! conquer
call fft_recursive(odd)
call fft_recursive(even)
! combine
do i=1,N/2
t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
x(i) = odd(i) + t
x(i+N/2) = odd(i) - t
end do
deallocate(odd)
deallocate(even)
end subroutine fft_recursive
subroutine initialize_filter
use space_grid
use array
use model
implicit none
integer :: i
integer :: N
integer :: center
N = (ize-izs+1)
center = int( ceiling( (N+1)*0.5 ) )
do i=izs,ize
fft_filter(i) = 1._dp- exp( -1._dp*(i-center)*(i-center)/ (2*fft_sigma*fft_sigma) ) ! or just a step function?
! write(*,'("(", F20.15, ",", F20.15, "i )")') fft_filter(i)
end do
end subroutine initialize_filter
subroutine suppress_high_freq(x)
use space_grid
use array
implicit none
integer :: i
real(dp), dimension(izs:ize), intent(inout) :: x
real(dp) :: invN
do i=izs,ize
fft_input(i) = x(i)
end do
call fft_recursive(fft_input) ! Fast Fourier Transform
do i=izs,ize ! Apply filter
fft_input(i) = fft_input(i)*fft_filter(i)
end do
call fft_recursive(fft_input) ! Invert the FFT
invN = 1._dp/(ize-izs+1._dp)
do i=izs,ize
x(i) = real(fft_input(ize-i+1))*invN
end do
end subroutine suppress_high_freq
end module fft
c4science · Help