source: trunk/00FlowSolve_PL/PROJECTS/NEW_SRC/filter4.f90 @ 4

subroutine filter
#define DEBUG_LEVEL 1
 use mpi_parameters
 use boundary_information
 use dependent_variables
 use counters_flags_etc
 implicit none

 integer                               :: sizeofreal,locnx
 include '../input/problem_size.h'
 include 'mpif.h'

! At designated time steps, apply fourth order compact spatial filters
if( mod(istep,i_filter) == 0  .and. istep .ne. 0 ) then

#if DEBUG_LEVEL >= 1
 if(myid==0) write(0,*) 'hello world from subroutine filter '
#endif

 if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then
   locnx = (nx-1)/nprocs
  else
   locnx = nx/nprocs
  endif

!!  create a derived data type
  call MPI_TYPE_EXTENT(MPI_REAL,sizeofreal,ierr)
  call MPI_TYPE_VECTOR(locnx,1,1,MPI_REAL,oneslice,ierr)
  call MPI_TYPE_HVECTOR(ny,1,nx*sizeofreal,oneslice,twoslice,ierr)
  call MPI_TYPE_HVECTOR(locnz,1,nx*ny*sizeofreal,twoslice,subslice,ierr)
  call MPI_TYPE_COMMIT(subslice,ierr)
  call MPI_BARRIER(comm,ierr)

  call filter4(U)
  call filter4(v)
  call filter4(W)
  call filter4(s1)
  call filter4(s2)

  call MPI_BARRIER(comm,ierr)

!! free data type
  call MPI_TYPE_FREE(oneslice,ierr)
  call MPI_TYPE_FREE(twoslice,ierr)
  call MPI_TYPE_FREE(subslice,ierr)

endif

end subroutine filter

subroutine filter4(U)
#define DEBUG_LEVEL 1
 use mpi_parameters
 use boundary_information
 use columnslab

 implicit none
 include 'mpif.h'

 include '../input/problem_size.h'
 integer                               :: n1,n2,n3,nzderiv,locnx
 real, dimension(:,:,:), allocatable   :: block_of_cols
 real, dimension(nx,ny,locnz)          :: U

 if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then
   locnx = (nx-1)/nprocs
  else
   locnx = nx/nprocs
  endif

! filter in x and y directions, overwriting input array in both cases
 call filter_x4(U,U,nx,ny,locnz)
 call filter_y4(U,U,nx,ny,locnz)

!filter in z direction :
!*******************************************************
! See compact_ddz_mpi for more details, but basically, each processor
! needs to ship sub-blocks of its own dataslice and recieve those
! necessary to construct a block with dimensions nx/numprocs * ny * nz.
! Filtering is done on these blocks and the results shipped back
! to the appropriate locations
!*******************************************************

 nzderiv=nz   ! don't use nzderiv < nz logic
!Both CASE1 and CASE2 require simple calls to compact_ddz
 if( numprocs == 1 .and. locnz == nz .and. nzderiv == nz ) then
  goto 999   ! CASE1
 elseif(numprocs > 1 .and. nzderiv == locnz .and. locnz < nz) then
  goto 999   ! CASE2
 elseif(numprocs > 1 .and. nzderiv == nz .and. nz > locnz) then
!CASE3 logic
 
  allocate ( block_of_cols( locnx,ny,numprocs*locnz) )

!!*******************************************************************************
!!*************   DO THE DATA EXCHANGES  ****************************************
!!*******************************************************************************

!!(1) Ship portions of my slice needed by other processes and
!!    receive data needed to fill in "block_of_cols"

   call slabs_to_columns(U,block_of_cols)

!!(2) Apply z4 filtering to "block_of_cols"

   n1=locnx
   n2=ny
   n3=locnz*numprocs
   call filter_z4(block_of_cols,block_of_cols,n1,n2,n3)

!!(3) Send filtered fields to other processes
!!    and receive their results.

   call columns_to_slabs(block_of_cols,U)
   
!!*******************************************************************************
!!*************    END DATA EXCHANGES    ****************************************
!!*******************************************************************************

  if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then
   U(nx,:,:)=U(1,:,:)  ! we didn't compute the values at extra bdry point location in x
  endif

  deallocate(block_of_cols)
  return
 endif

999 continue   ! CASE1 & CASE2
 n1=nx
 n2=ny
 n3=nzderiv
 call filter_z4(U,U,n1,n2,n3)

end subroutine filter4

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine filter_x4(data,deriv,nx,ny,nz)
! routine to perform 4th order filtering wrt x (1st dimension)
! of a 3d field stored in data(1:nx,1:ny,1:nz).
!
! inputs
!        data       nx,ny,nz array of input values
!        nx,ny,nz   dimensions of field to be differentiated
!
! output
!        deriv      nx,ny,nz array containing the 4th order
!                   filtered (in x) field
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 use grid_info
 use boundary_information
use user_parameters,        only: LAPACK_FLAG
 implicit none

!!! double precision, dimension(:,:),allocatable    :: work
 real, dimension(:,:),allocatable    :: work
 integer                             :: i,j,k,nx,ny,nz,ierr
 real                                :: data(nx,ny,nz),deriv(nx,ny,nz)
!!! double precision                    :: a,b,c
 real                    :: a,b,c

 real,allocatable,dimension(:)       :: kx,SF,in,out,tmp
 real                                :: pi,kx_max,xnorm,kstar,dk
 integer*8                           :: plan_f,plan_i     !size must match pointer size, 8 for DEC alphas
 integer                             :: N
 include 'fftw_kw.inc'               !  predefined FFTW constants


 if( boundary_type2(1,1)=='periodic' ) then
  N=nx-1                 ! nx includes both boundary points, N only 1
  allocate(kx(N+1),SF(N))
  allocate( in(N),out(N),tmp(N) )

! X preliminaries, setup fftw "plan"
  call rfftw_f77_create_plan(plan_f,N,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE)
  call rfftw_f77_create_plan(plan_i,N,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE)


! create the x array, align with the data arrays
  pi=4.*atan(1.)
  xnorm = 1./float(N)    ! only normalize in wavenumber space here
  do i=1,N/2+1
   kx(i)=2.*pi*(i-1.)    ! FFTW real to half-complex storage
  enddo
  do i=2,N/2
   kx(N-i+2)=kx(i)       ! FFTW real to half-complex storage
  enddo
  kx(N+1)=0.             ! extra location, convenient for alignment of arrays

! define smoothing filter
  kx_max=kx(N/2+1)
  kstar=(8./9.)*kx_max 
  do i=1,N
       !if( kx(i) >= (8./9.)*kx_max ) then
       ! SF(i)=0.
       !else
       ! SF(i)=1.0
       !endif
   SF(i)=0.5*( -tanh( 8.*(kx(i)-kstar)/(kx_max) ) + 1. )
  enddo

  do k=1,nz
   do j=1,ny

! x transforms; filter & inverse transform
    in(1:N)=data(1:N,j,k)
    call rfftw_f77_one(plan_f,in,out)
    tmp(1:N)=SF(1:N)*out(1:N)*xnorm
    call rfftw_f77_one(plan_i,tmp,out)
    deriv(1:N,j,k)=out
    deriv(nx,j,k)=deriv(1,j,k)

   enddo
  enddo

  call rfftw_f77_destroy_plan(plan_f)
  call rfftw_f77_destroy_plan(plan_i)
  deallocate( kx,SF,in,out,tmp )
  return
 endif

!               GENERAL, NONPERIODIC LOGIC
allocate ( work(nx,ny) )
!Set filter coefficients:
 a=(5.+6.*alpha)/8.
 b=(1.+2.*alpha)/2.
 c=-(1.-2.*alpha)/8.

!Loop through each horizontal plane
 do k=1,nz

!Create rhs vectors, store in work(:,:,1)
  do i=3,nx-2
   do j=1,ny
    work(i,j)=a*data(i,j,k) + b*( data(i+1,j,k)+data(i-1,j,k) )/2.   &
               + c*( data(i+2,j,k)+data(i-2,j,k) )/2.
   enddo
  enddo

  do j=1,ny
   work(1,j)= (15./16.)*data(1,j,k) + (1./16.)*( 4*data(2,j,k)         &
               -6.*data(3,j,k) + 4.*data(4,j,k) - data(5,j,k) )
   work(2,j)= (3./4.)*data(2,j,k) + (1./16.)*( data(1,j,k)             &
               +6.*data(3,j,k) - 4.*data(4,j,k) + data(5,j,k) )

   work(nx,j)=(15./16.)*data(nx,j,k) + (1./16.)*( 4*data(nx-1,j,k)     &
               -6.*data(nx-2,j,k) + 4.*data(nx-3,j,k) - data(nx-4,j,k))
   work(nx-1,j)= (3./4.)*data(nx-1,j,k) + (1./16)*( data(nx,j,k)      &
               +6.*data(nx-2,j,k) - 4.*data(nx-3,j,k) + data(nx-4,j,k))
  enddo

! compute the filtered fields for plane k, results overwritten onto work(:,:)
 if( LAPACK_FLAG == 'double' ) then
!!!  call DGBTRS('N',nx,1,1,ny,Grid(0)%filt_x,4,Grid(0)%piv_fx,work(1,1),nx,ierr)
      write(0,*) 'LAPACK_FLAG = double instead of single'
      stop
 elseif( LAPACK_FLAG == 'single' ) then
  call SGBTRS('N',nx,1,1,ny,Grid(0)%filt_x,4,Grid(0)%piv_fx,work(1,1),nx,ierr)
 endif
  
  deriv(:,:,k)=work(:,:)     !store the results in deriv(:,:,k)

 enddo ! next k plane

 deallocate ( work )
end subroutine filter_x4

subroutine filter_y4(data,deriv,nx,ny,nz)
! routine to perform 4th order filtering wrt y (2nd dimension)
! of a 3d field stored in data(1:nx,1:ny,1:nz).
!
! inputs
!        data       nx,ny,nz array of input values
!        nx,ny,nz   dimensions of field to be differentiated

! output
!        deriv      nx,ny,nz array containing the 4th order
!                   accurate filtered (in y ) field
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 use grid_info
 use boundary_information
use user_parameters,        only: LAPACK_FLAG
 implicit none

!!! double precision, dimension(:,:),allocatable    :: work
 real, dimension(:,:),allocatable    :: work
 integer                             :: i,j,k,nx,ny,nz,ierr
 real                                :: data(nx,ny,nz),deriv(nx,ny,nz)
!!! double precision                    :: a,b,c
 real                    :: a,b,c

 real,allocatable,dimension(:)       :: ky,SF,in,out,tmp
 real                                :: pi,ky_max,xnorm,kstar,dk
 integer*8                           :: plan_f,plan_i     !size must match pointer size, 8 for DEC alphas
 integer                             :: N
 include 'fftw_kw.inc'               !  predefined FFTW constants


 if( boundary_type4(1,1)=='periodic' ) then
  N=ny-1                 ! ny includes both boundary points, N only 1
  allocate( in(N),out(N),tmp(N) )
  allocate(ky(N+1),SF(N))

! Y preliminaries, setup fftw "plan"
  call rfftw_f77_create_plan(plan_f,N,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE)
  call rfftw_f77_create_plan(plan_i,N,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE)

! create the y array, align with the data arrays
  pi=4.*atan(1.)
  xnorm = 1./float(N)    ! only normalize in wavenumber space here
  do j=1,N/2+1
   ky(j)=2.*pi*(j-1.)    ! FFTW real to half-complex storage
  enddo
  do j=2,N/2
   ky(N-j+2)=ky(j)       ! FFTW real to half-complex storage
  enddo
  ky(N+1)=0.             ! extra location, convenient for alignment of arrays

! define smoothing filter
  ky_max=ky(N/2+1)
  kstar=(8./9.)*ky_max

  if( N <= 16 ) then
   do j=1,N
     SF(j)=1.0    ! no filtering
   enddo
  elseif( N > 16 ) then
   do j=1,N
    SF(j)=0.5*( -tanh( 8.*(ky(j)-kstar)/(ky_max) ) + 1. )
   enddo
  endif

  do k=1,nz
   do i=2,nx-1

! y transforms; filter & inverse transform
    in(1:N)=data(i,1:N,k)
    call rfftw_f77_one(plan_f,in,out)
    tmp(1:N)=SF(1:N)*out(1:N)*xnorm
    call rfftw_f77_one(plan_i,tmp,out)
    deriv(i,1:N,k)=out
    deriv(i,ny,k)=deriv(i,1,k)

   enddo
  enddo
  deriv(1,:,:)=data(1,:,:)
  deriv(nx,:,:)=data(nx,:,:)

  call rfftw_f77_destroy_plan(plan_f)
  call rfftw_f77_destroy_plan(plan_i)
  deallocate( ky,SF,in,out,tmp )
  return
 endif

!               GENERAL, NONPERIODIC LOGIC

 allocate ( work(ny,nx) )
 a=(5.+6.*alpha)/8.
 b=(1.+2.*alpha)/2.
 c=-(1.-2.*alpha)/8.

!Loop through each horizontal plane
 do k=1,nz
! create rhs vectors, store in work(:,:,1)
  do j=3,ny-2
   do i=1,nx
    work(j,i)=a*data(i,j,k) + b*( data(i,j+1,k)+data(i,j-1,k) )/2. &
               + c*( data(i,j+2,k)+data(i,j-2,k) )/2.
   enddo
  enddo

  do i=1,nx
   work(1,i)= (15./16.)*data(i,1,k) + (1./16.)*( 4.*data(i,2,k)          &
              -6.*data(i,3,k) + 4.*data(i,4,k) - data(i,5,k) )
   work(2,i)= (3./4.)*data(i,2,k) + (1./16.)*( data(i,1,k)              &
              +6.*data(i,3,k) - 4.*data(i,4,k) + data(i,5,k) )

   work(ny,i)=(15./16.)*data(i,ny,k) + (1./16.)*( 4.*data(i,ny-1,k)      &
               -6.*data(i,ny-2,k) + 4.*data(i,ny-3,k) - data(i,ny-4,k))
   work(ny-1,i)=(3./4.)*data(i,ny-1,k) + (1./16)*( data(i,ny,k)        &
                 +6.*data(i,ny-2,k) - 4.*data(i,ny-3,k) + data(i,ny-4,k))
  enddo

! compute the filtered fields for plane k, overwriting results into work(:,:)
 if( LAPACK_FLAG == 'double' ) then
!!!  call DGBTRS('N',ny,1,1,nx,Grid(0)%filt_y,4,Grid(0)%piv_fy,work(1,1),ny,ierr)
      write(0,*) 'LAPACK_FLAG = double instead of single'
      stop
 elseif( LAPACK_FLAG == 'single' ) then
  call SGBTRS('N',ny,1,1,nx,Grid(0)%filt_y,4,Grid(0)%piv_fy,work(1,1),ny,ierr)
 endif

!     store the results in deriv(:,:,k)
  do i=1,nx
   do j=1,ny
    deriv(i,j,k)=work(j,i)
   enddo
  enddo

 enddo   ! next k plane

 deallocate ( work )
end subroutine filter_y4

subroutine filter_z4(data,deriv,nx,ny,nz)
!
!       routine to perform 4th order filtering wrt z (3rd dimension)
!       of a 3d field stored in data(1:nx,1:ny,1:nz).
!
! inputs
!        data       nx,ny,nz array of input values
!        nx,ny,nz   dimensions of field to be filtered
!
! output
!        deriv      nx,ny,nz array containing the filtered field
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 use grid_info
 use user_parameters, only: LAPACK_FLAG
 implicit none

!!! double precision, dimension(:,:),allocatable    :: work
 real, dimension(:,:),allocatable    :: work
 integer                             :: i,j,k,nx,ny,nz,ierr
 real                                :: data(nx,ny,nz),deriv(nx,ny,nz)
!!! double precision                    :: a,b,c
 real                    :: a,b,c
 allocate ( work(nz,nx) )
      
 a=(5.+6.*alpha)/8.
 b=(1.+2.*alpha)/2.
 c=-(1.-2.*alpha)/8.

!Loop through each xz vertical plane
 do j=1,ny
! create rhs vectors, store in work(:,:,1)

  do k=3,nz-2
   do i=1,nx
    work(k,i)=a*data(i,j,k) + b*( data(i,j,k+1)+data(i,j,k-1) )/2.   &
               + c*( data(i,j,k+2)+data(i,j,k-2) )/2.
   enddo
  enddo

  do i=1,nx
   work(1,i)= (15./16.)*data(i,j,1) + (1./16.)*( 4.*data(i,j,2)             &
              -6.*data(i,j,3) + 4.*data(i,j,4) - data(i,j,5) )
   work(2,i)= (3./4.)*data(i,j,2) + (1./16.)*( data(i,j,1)                 &
              +6.*data(i,j,3) - 4.*data(i,j,4) + data(i,j,5) )

   work(nz,i)=(15./16.)*data(i,j,nz) + (1./16.)*( 4.*data(i,j,nz-1)         &
               -6.*data(i,j,nz-2) + 4.*data(i,j,nz-3) - data(i,j,nz-4))
   work(nz-1,i)=(3./4.)*data(i,j,nz-1) + (1./16.)*( data(i,j,nz)           &
                 +6.*data(i,j,nz-2) - 4.*data(i,j,nz-3) + data(i,j,nz-4))
  enddo

!Filter the fields wrt z for plane j, overwriting results into work(:,:)
 if( LAPACK_FLAG == 'double' ) then
!!!  call DGBTRS('N',nz,1,1,nx,Grid(0)%filt_z,4,Grid(0)%piv_fz,work(1,1),nz,ierr)
      write(0,*) 'LAPACK_FLAG = double instead of single'
      stop
 elseif( LAPACK_FLAG == 'single' ) then
  call SGBTRS('N',nz,1,1,nx,Grid(0)%filt_z,4,Grid(0)%piv_fz,work(1,1),nz,ierr)
 endif

!Store the results in deriv(:,j,:)
  do k=1,nz
   do i=1,nx
    deriv(i,j,k)=work(k,i)
   enddo
  enddo

 enddo    ! go on to the next plane

 deallocate ( work )
end subroutine filter_z4