poisson
module decomp_2d_poisson
use decomp_2d
use decomp_2d_fft
use param
use variables
implicit none
private ! Make everything private unless declared public
! real(mytype), private, parameter :: PI = 3.14159265358979323846_mytype
! This 0.0_mytype!!!!!!!!!!!!!1111
#ifdef DOUBLE_PREC
real(mytype), parameter :: epsilon = 1.e-16
#else
real(mytype), parameter :: epsilon = 1.e-8
#endif
! boundary conditions
integer, save :: bcx, bcy, bcz
! decomposition object for physical space
TYPE(DECOMP_INFO), save :: ph
! decomposition object for spectral space
TYPE(DECOMP_INFO), save :: sp
! store sine/cosine factors
real(mytype), save, allocatable, dimension(:) :: az,bz
real(mytype), save, allocatable, dimension(:) :: ay,by
real(mytype), save, allocatable, dimension(:) :: ax,bx
! wave numbers
complex(mytype), save, allocatable, dimension(:,:,:) :: kxyz
!wave numbers for stretching in a pentadiagonal matrice
complex(mytype), save, allocatable, dimension(:,:,:,:) :: a,a2,a3
! work arrays,
! naming convention: cw (complex); rw (real);
! 1 = X-pencil; 2 = Y-pencil; 3 = Z-pencil
real(mytype), allocatable, dimension(:,:,:) :: rw1,rw1b,rw2,rw2b,rw3
complex(mytype), allocatable, dimension(:,:,:) :: cw1,cw1b,cw2,cw22,cw2b,cw2c
! underlying FFT library only needs to be initialised once
logical, save :: fft_initialised = .false.
public :: decomp_2d_poisson_stg, decomp_2d_poisson_init, &
decomp_2d_poisson_finalize
! For staggered mesh where main variables are defined in the centre of
! control volumes while boundary conditions are defined on interfaces
interface decomp_2d_poisson_stg
module procedure poisson
end interface
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialise Poisson solver for given boundary conditions
! given
! given
! just for init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine decomp_2d_poisson_init(bcx1, bcy1, bcz1)
implicit none
integer, intent(IN) :: bcx1, bcy1, bcz1
integer :: nx, ny, nz, i
bcx = bcx1
bcy = bcy1
bcz = bcz1
nx = nx_global
ny = ny_global
nz = nz_global
! pressure-grid having 1 fewer point for non-periodic directions
if (bcx==1) nx=nx-1
if (bcy==1) ny=ny-1
if (bcz==1) nz=nz-1
allocate(ax(nx),bx(nx))
allocate(ay(ny),by(ny))
allocate(az(nz),bz(nz))
call abxyz(ax,ay,az,bx,by,bz,nx,ny,nz,bcx,bcy,bcz)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Initialise the ax bx
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Initialise the ay by
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Initialise the az bz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
call decomp_info_init(nx, ny, nz, ph)
call decomp_info_init(nx, ny, nz/2+1, sp)
! allocate work space
if (bcx==0 .and. bcy==0 .and. bcz==0) then
allocate(cw1(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(kxyz(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(a(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a2(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a3(sp%yst(1):sp%yen(1),ny,sp%yst(3):sp%yen(3),5))
! (000) ^c ^kxyz ^a
else if (bcx==1 .and. bcy==0 .and. bcz==0) then
allocate(cw1(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(cw1b(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(rw1(ph%xst(1):ph%xen(1),ph%xst(2):ph%xen(2), &
ph%xst(3):ph%xen(3)))
allocate(rw1b(ph%xst(1):ph%xen(1),ph%xst(2):ph%xen(2), &
ph%xst(3):ph%xen(3)))
allocate(rw2(ph%yst(1):ph%yen(1),ph%yst(2):ph%yen(2), &
ph%yst(3):ph%yen(3)))
allocate(kxyz(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(a(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a2(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a3(sp%yst(1):sp%yen(1),ny,sp%yst(3):sp%yen(3),5))
!(100) consist :: cw1 cw1b rw1 rw1b rw2 ^a
else if (bcx==0 .and. bcy==1 .and. bcz==0) then
allocate(rw2(ph%yst(1):ph%yen(1),ph%yst(2):ph%yen(2), &
ph%yst(3):ph%yen(3)))
allocate(rw2b(ph%yst(1):ph%yen(1),ph%yst(2):ph%yen(2), &
ph%yst(3):ph%yen(3)))
allocate(cw1(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(cw2(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(cw22(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(cw2b(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(cw2c(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(kxyz(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(a(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a2(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a3(sp%yst(1):sp%yen(1),ny,sp%yst(3):sp%yen(3),5))
!(010) rw2 rw2b cw2 cw22 cw2b cw2c kxyz ^a
else if (bcx==1 .and. bcy==1) then
allocate(cw1(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(cw1b(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(cw2(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(cw22(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(cw2b(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(cw2c(sp%yst(1):sp%yen(1),sp%yst(2):sp%yen(2), &
sp%yst(3):sp%yen(3)))
allocate(rw1(ph%xst(1):ph%xen(1),ph%xst(2):ph%xen(2), &
ph%xst(3):ph%xen(3)))
allocate(rw1b(ph%xst(1):ph%xen(1),ph%xst(2):ph%xen(2), &
ph%xst(3):ph%xen(3)))
allocate(rw2(ph%yst(1):ph%yen(1),ph%yst(2):ph%yen(2), &
ph%yst(3):ph%yen(3)))
allocate(rw2b(ph%yst(1):ph%yen(1),ph%yst(2):ph%yen(2), &
ph%yst(3):ph%yen(3)))
!(11*) cw1 cw1b cw2 cw2b cw2c rw1 rw1b rw2b
if (bcz==1) then
allocate(rw3(ph%zsz(1),ph%zsz(2),ph%zsz(3)))
end if
allocate(kxyz(sp%xst(1):sp%xen(1),sp%xst(2):sp%xen(2), &
sp%xst(3):sp%xen(3)))
allocate(a(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a2(sp%yst(1):sp%yen(1),ny/2,sp%yst(3):sp%yen(3),5))
allocate(a3(sp%yst(1):sp%yen(1),nym,sp%yst(3):sp%yen(3),5))
!(111) rw3 kxyz ^a
end if
call waves()
return
end subroutine decomp_2d_poisson_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Release memory used by Poisson solver
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine decomp_2d_poisson_finalize
implicit none
deallocate(ax,bx,ay,by,az,bz)
call decomp_info_finalize(ph)
call decomp_info_finalize(sp)
call decomp_2d_fft_finalize
fft_initialised = .false.
deallocate(kxyz)
if (bcx==0 .and. bcy==0 .and. bcz==0) then
deallocate(cw1)
! (000) cw1
else if (bcx==1 .and. bcy==0 .and. bcz==0) then
deallocate(cw1,cw1b,rw1,rw1b,rw2)
! (100) cw1 cw1b rw1 rw1b rw2
else if (bcx==0 .and. bcy==1 .and. bcz==0) then
deallocate(cw1,cw2,cw2b,rw2,rw2b)
! (010) cw1 cw2 cw2b rw2 rw2b
else if (bcx==1 .and. bcy==1) then
deallocate(cw1,cw1b,cw2,cw2b,rw1,rw1b,rw2,rw2b)
! (11*) cw1,cw1b,cw2,cw2b,rw1,rw1b,rw2,rw2b
if (bcz==1) then
deallocate(rw3)
!(111) rw3
end if
end if
return
end subroutine decomp_2d_poisson_finalize
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Top level wrapper
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine poisson(rhs, bcx, bcy, bcz)
implicit none
real(mytype), dimension(:,:,:), intent(INOUT) :: rhs
integer, intent(IN) :: bcx, bcy, bcz ! boundary conditions
integer :: i
if (bcx==0 .and. bcy==0 .and. bcz==0) then
call poisson_000(rhs)
else if (bcx==1 .and. bcy==0 .and. bcz==0) then
call poisson_100(rhs)
else if (bcx==0 .and. bcy==1 .and. bcz==0) then
call poisson_010(rhs)
else if (bcx==1 .and. bcy==1) then ! 110 & 111
call poisson_11x(rhs, bcz)
else
stop 'boundary condition not supported'
end if
return
end subroutine poisson
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solving 3D Poisson equation with periodic B.C in all 3 dimensions
! 3
! 3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine poisson_000(rhs)
use derivX
use derivY
use derivZ
! right-hand-side of Poisson as input
! solution of Poisson as output
real(mytype), dimension(:,:,:), intent(INOUT) :: rhs
integer, dimension(3) :: fft_start, fft_end, fft_size
complex(mytype) :: xyzk
complex(mytype) :: ytt,xtt,ztt,yt1,xt1,yt2,xt2
complex(mytype) :: xtt1,ytt1,ztt1,zt1,zt2
real(mytype) :: tmp1, tmp2,x ,y, z
integer :: nx,ny,nz, i,j,k
nx = nx_global
ny = ny_global
nz = nz_global
if (.not. fft_initialised) then
call decomp_2d_fft_init(PHYSICAL_IN_Z)
fft_initialised = .true.
end if
! compute r2c transform r2c:real to complex by zhaoliang
! cw1 is 3dim data structrue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! only fft to transform the real data to complex data while (100) twice
!!!!! one is tranform and then fft (you can see poisson100
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call decomp_2d_fft_3d(rhs,cw1)
! normalisation
! why this below is called normalisation ?????? just divide 3dim grid number
cw1 = cw1 / real(nx, kind=mytype) /real(ny, kind=mytype) &
/ real(nz, kind=mytype)
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
! post-processing in spectral space
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!post-processing in 3-direction can be
!!!!!!!!!!!!!!!!! set in the same cycle,because they are
!!!!!!!!!!!!!!!!!! similar
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! POST PROCESSING IN Z
tmp1 = real(cw1(i,j,k), kind=mytype) ! get the real number of cw1
tmp2 = aimag(cw1(i,j,k)) ! get the imagenumber of cw1
cw1(i,j,k) = cmplx(tmp1*bz(k)+tmp2*az(k), &
tmp2*bz(k)-tmp1*az(k), kind=mytype) ! modify the cw1 in the spectral space Z
! bz az
! by ay
! bx ax
! POST PROCESSING IN Y
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*by(j)+tmp2*ay(j), &
tmp2*by(j)-tmp1*ay(j), kind=mytype)
if (j.gt.(ny/2+1)) cw1(i,j,k)=-cw1(i,j,k) ! why should be axisymmetry!
! POST PROCESSING IN X
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bx(i)+tmp2*ax(i), &
tmp2*bx(i)-tmp1*ax(i), kind=mytype)
if (i.gt.(nx/2+1)) cw1(i,j,k)=-cw1(i,j,k)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!! Solve Poisson
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
tmp1=real(kxyz(i,j,k), kind=mytype)
tmp2=aimag(kxyz(i,j,k))
! CANNOT DO A DIVISION BY ZERO
! Yes ! division by zero is impossible!!!-----------------------------------<
if ((tmp1.lt.epsilon).or.(tmp2.lt.epsilon)) then !epsilon? what does it mean?
cw1(i,j,k)=0._mytype
! print *,'DIV 0',i,j,k,epsilon
else
cw1(i,j,k)=cmplx( real(cw1(i,j,k), kind=mytype) / (-tmp1), &
aimag(cw1(i,j,k))/(-tmp2), kind=mytype)
end if
!Print result in spectal space after Poisson
! if (abs(out(i,j,k)) > 1.0e-4) then
! write(*,*) 'AFTER',i,j,k,out(i,j,k),xyzk
! end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! post-processing backward
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! POST PROCESSING IN Z
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)-tmp2*az(k), &
-tmp2*bz(k)-tmp1*az(k), kind=mytype)
! bz az
! by ay
! bx ax
! POST PROCESSING IN Y
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*by(j)+tmp2*ay(j), &
tmp2*by(j)-tmp1*ay(j), kind=mytype)
if (j.gt.(ny/2+1)) cw1(i,j,k)=-cw1(i,j,k)
! POST PROCESSING IN X
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bx(i)+tmp2*ax(i), &
-tmp2*bx(i)+tmp1*ax(i), kind=mytype)
if (i.gt.(nx/2+1)) cw1(i,j,k)=-cw1(i,j,k)
end do
end do
end do
! compute c2r transform
call decomp_2d_fft_3d(cw1,rhs) ! from complex to real!
! call decomp_2d_fft_finalize
return
end subroutine poisson_000
subroutine poisson_100(rhs)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! get the rhs and then modify rhs ,at last return the rhs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
real(mytype), dimension(:,:,:), intent(INOUT) :: rhs
complex(mytype) :: xyzk
real(mytype) :: tmp1, tmp2, tmp3, tmp4
real(mytype) :: xx1,xx2,xx3,xx4,xx5,xx6,xx7,xx8
integer :: nx,ny,nz, i,j,k, itmp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! 8 transform from x -->y -->z then z--> y --->x at the beginning
!!!!!!!!!!!!!!!!!!from x -->y -->z then z--> y --->x at the endding
!!!!! (100 add 8 more transform than 000) while 000 is 0 transform
!!!!! so rhs is in the z pencil ,no there are not z pencil at all!
!!!!! the program is 2-dim pm*pn decomp ! so only xpencil and y pencil
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!twice real to complex
!!!!!!!!!!! 1: transform rhs-rw2-rw1
!!!!!!!!!!! 2: fft rhs-----rw1
!!!!!!!!!!!!!!twice complex to real
!!!!!!!!!!! 1: transfrom rw1-rw2--rhs
!!!!!!!!!!! 2: fft rw1------rhs
!!!!!!!!!!!!! post-processing operations is under the complex condition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!! input real rhs(:.:,:) in the z pencil
!!!!!!!!!!!! 2+2 transform |destination: modify the rhs
!!!!!!!!!!!! 1 normalisation
!!!!!!!!!!!! 1 FFT forward
!!!!!!!!!!!! 3 post-processing z-->y--->x direction not pencil
!!!!!!!!!!!! 1 poisson solver
!!!!!!!!!!!! 3 post-processing backward x--->y---->z direction
!!!!!!!!!!!! 1 FFT backward
!!!!!!!!!!!! 2+2 transform |destination: modify the rhs
!!!!!!!!!!!! output real rhs(:,:,:) in the z pencil
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
100 format(1x,a8,3I4,2F12.6)
nx = nx_global - 1
ny = ny_global
nz = nz_global
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
! rhs is in Z-pencil(wrong!There isn't Z pencil)(wrong
!again, there is Z pencil))but requires global operations in X
! why should do the transpose!!!!!!!!!!!!!!
! why in the poisson_000 is not needed!???????????????
! two steps to change from z to x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
call transpose_z_to_y(rhs,rw2,ph) ! z pencil in y pencil but in the physical space
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!rw2 is the complex 3-dim data the same as rw1
call transpose_y_to_x(rw2,rw1,ph)
do k=ph%xst(3),ph%xen(3)
do j=ph%xst(2),ph%xen(2)
do i=1,nx/2
!!! rw1 rw1b 1 means x pencil
!!! rw2 rw2b 2 means y pencil
rw1b(i,j,k)=rw1(2*(i-1)+1,j,k) ! the odd terms is for the half before
enddo
do i=nx/2+1,nx
rw1b(i,j,k)=rw1(2*nx-2*i+2,j,k) ! the reverse terms is for the half after
enddo
enddo
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
!!!!3 layers data
!!!! rw1b is the top lay data xpencil physical space
!!!! rw2 is in the middle lay data y pencil physical space
!!!! rw1 rhs is in the outer lay data z pencil spectral space
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
call transpose_x_to_y(rw1b,rw2,ph) ! let x pencil rw1b to y pencil
call transpose_y_to_z(rw2,rhs,ph) ! let y pencil rw2 to z pencil's rhs
!!!!!!!!!!!!first time get what we want
if (.not. fft_initialised) then
call decomp_2d_fft_init(PHYSICAL_IN_Z,nx,ny,nz)
fft_initialised = .true.
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! compute r2c transform
! why should force the before pre-posting ------------------------------<<<<
!!!!!!!!!the most key operation :::: the fft operator: decomp_2d_fft_3d
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call decomp_2d_fft_3d(rhs,cw1) ! from real to complex! and begin poission solver
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! in the z pencil we did the fft transform
!!!!!!!!!!!!!!!!!! so cw1 is 3d complex data ,yeah (:,:,:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! normalisation
cw1 = cw1 / real(nx, kind=mytype) /real(ny, kind=mytype) &
/ real(nz, kind=mytype)
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'START',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! post-processing in spectral space
! POST PROCESSING IN Z
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)+tmp2*az(k), &
tmp2*bz(k)-tmp1*az(k), kind=mytype)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'after z',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! POST PROCESSING IN Y
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*by(j)+tmp2*ay(j), &
tmp2*by(j)-tmp1*ay(j), kind=mytype)
if (j.gt.(ny/2+1)) cw1(i,j,k)=-cw1(i,j,k)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'after y',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!because (100)
!!!!!!!!!!!!!!!!!!!!!so the post processing in x is different
!!!!!!!!!!!!!!!!!!!!!from the y and z
! POST PROCESSING IN X
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
cw1b(1,j,k)=cw1(1,j,k)
do i = 2,nx
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
tmp3 = real(cw1(nx-i+2,j,k), kind=mytype) !!! the 2~nx
tmp4 = aimag(cw1(nx-i+2,j,k)) !!!!!!!!!!!! 2~nx
xx1=tmp1*bx(i)/2._mytype
xx2=tmp1*ax(i)/2._mytype
xx3=tmp2*bx(i)/2._mytype
xx4=tmp2*ax(i)/2._mytype
xx5=tmp3*bx(i)/2._mytype
xx6=tmp3*ax(i)/2._mytype
xx7=tmp4*bx(i)/2._mytype
xx8=tmp4*ax(i)/2._mytype
cw1b(i,j,k) = cmplx(xx1+xx4+xx5-xx8,-xx2+xx3+xx6+xx7, &
kind=mytype)
end do
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! While you set the -Debug in the compile period
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1b(i,j,k)) > 1.0e-4) then
write(*,100) 'after x',i,j,k,cw1b(i,j,k)
end if
end do
end do
end do
#endif
!!!!!!!!!!!!!!!!!!!!!from now on, the cw1b in x,y,z direction(not pencil)
!!!!!!!!!!!!!!!!!!!!is calculated ,so now you can calculate cw1b in the
!!!!!!!!!!!!!!!!!!!!!!! poisson condition ,yes poisson solver
! Solve Poisson
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
!tmp1=real(zk2(k)+yk2(j)+xk2(i), kind=mytype)
!tmp2=aimag(zk2(k)+yk2(j)+xk2(i))
tmp1=real(kxyz(i,j,k), kind=mytype)
tmp2=aimag(kxyz(i,j,k))
!xyzk=cmplx(tmp1,tmp2, kind=mytype)
! CANNOT DO A DIVISION BY ZERO
! yes
if ((abs(tmp1).lt.epsilon).and.(abs(tmp2).lt.epsilon)) then
cw1b(i,j,k)=cmplx(0._mytype,0._mytype, kind=mytype)
end if
if ((abs(tmp1).lt.epsilon).and.(abs(tmp2).ge.epsilon)) then
cw1b(i,j,k)=cmplx(0._mytype, &
aimag(cw1b(i,j,k))/(-tmp2), kind=mytype)
end if
if ((abs(tmp1).ge.epsilon).and.(abs(tmp2).lt.epsilon)) then
cw1b(i,j,k)=cmplx( real(cw1b(i,j,k), kind=mytype) &
/(-tmp1), 0._mytype, kind=mytype)
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!the most dayly processing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if ((abs(tmp1).ge.epsilon).and.(abs(tmp2).ge.epsilon)) then
cw1b(i,j,k)=cmplx( real(cw1b(i,j,k), kind=mytype) &
/(-tmp1), &
aimag(cw1b(i,j,k))/(-tmp2), kind=mytype)
end if
#ifdef DEBUG
if (abs(cw1b(i,j,k)) > 1.0e-4) &
write(*,100) 'AFTER',i,j,k,cw1b(i,j,k)
#endif
end do
end do
end do
! post-processing backward
!!!!!!!!!!!!!!!!!!!!!! You know cw1b is the cw1's another form in the spectral
!!!!!!!!!!!!!!!!!!!!!! space ,so you need to change the spetral space to physical
!!!!!!!!!!!!!!!!!!!!!!! space again ,yes get the cw1
!!!!!!!!!!!!!!!!!!!!!!!! from now on ,you will found that the calculate of poisson
!!!!!!!!!!!!!!!!!!!!!!!! equation is in the spectral space ,from the fft,then 3 post
!!!!!!!!!!!!!!!!!!!!!!!!!! processing in x,y,z 3 directions, then calculate the poisson
!!!!!!!!!!!!!!!!!!!!!!!!!! equation. OK,then you start next!
! POST PROCESSING IN X
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
cw1(1,j,k)=cw1b(1,j,k)
do i = 2,nx
tmp1 = real(cw1b(i,j,k), kind=mytype)
tmp2 = aimag(cw1b(i,j,k))
tmp3 = real(cw1b(nx-i+2,j,k), kind=mytype)
tmp4 = aimag(cw1b(nx-i+2,j,k))
xx1=tmp1*bx(i)
xx2=tmp1*ax(i)
xx3=tmp2*bx(i)
xx4=tmp2*ax(i)
xx5=tmp3*bx(i)
xx6=tmp3*ax(i)
xx7=tmp4*bx(i)
xx8=tmp4*ax(i)
cw1(i,j,k) = cmplx(xx1-xx4+xx6+xx7,-(-xx2-xx3+xx5-xx8), &
kind=mytype)
!!!!!the buterfly algorithm!!!!!!!!!!!!!!!!!!!!!!!
end do
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!from now on x direction ,the spectral data has been transform to
!!!!!!!!!!!!!!!!!!!!!!!! physical space, but still the complex
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'AFTER X',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! POST PROCESSING IN Y
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*by(j)-tmp2*ay(j), &
tmp2*by(j)+tmp1*ay(j), kind=mytype)
if (j.gt.(ny/2+1)) cw1(i,j,k)=-cw1(i,j,k)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'AFTER Y',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
!!!!!!!!!!!!!!!!!!!!!! by ay
!!!!!!!!!!!!!!!!!!!!!! bz az
! POST PROCESSING IN Z
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)-tmp2*az(k), &
tmp2*bz(k)+tmp1*az(k), kind=mytype)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'END',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! compute c2r transform
call decomp_2d_fft_3d(cw1,rhs)
! rhs is in Z-pencil but requires global operations in X
call transpose_z_to_y(rhs,rw2,ph)
call transpose_y_to_x(rw2,rw1,ph)
do k=ph%xst(3),ph%xen(3)
do j=ph%xst(2),ph%xen(2)
do i=1,nx/2
rw1b(2*i-1,j,k)=rw1(i,j,k) ! the odd terms is setted!
enddo
do i=1,nx/2
rw1b(2*i,j,k)=rw1(nx-i+1,j,k) ! the even terms is setted
enddo
enddo
end do
call transpose_x_to_y(rw1b,rw2,ph)
call transpose_y_to_z(rw2,rhs,ph)
!!!!!!!!!!!!!the finally outcome rhs
! call decomp_2d_fft_finalize
return
end subroutine poisson_100
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Neumann!!!!!!!!!!!!!!!!!!!!!!!!!!!111
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solving 3D Poisson equation: Neumann in Y; periodic in X & Z
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine poisson_010(rhs)
implicit none
real(mytype), dimension(:,:,:), intent(INOUT) :: rhs
complex(mytype) :: xyzk
real(mytype) :: tmp1, tmp2, tmp3, tmp4
real(mytype) :: xx1,xx2,xx3,xx4,xx5,xx6,xx7,xx8
integer :: nx,ny,nz, i,j,k
100 format(1x,a8,3I4,2F12.6)
nx = nx_global
ny = ny_global - 1 ! (010) so have fewer1 ny!
nz = nz_global
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!so if (11*) you need to ?what pencil??????
!!!!??
! rhs is in Z-pencil but requires global operations in Y because(010)
!!!!!!!!!!!!!!twice real to complex
!!!!!!!!!!! 1: transform rhs-rw2-rw1
!!!!!!!!!!! 2: fft rhs-----rw1
!!!!!!!!!!!!!!twice complex to real
!!!!!!!!!!! 1: transfrom rw1-rw2--rhs
!!!!!!!!!!! 2: fft rw1------rhs
!!!!!!!!!!!!! post-processing operations is under the complex condition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!! input real rhs(:.:,:) in the z pencil
!!!!!!!!!!!! 1+1 transform z->y y->z |destination: modify the rhs because z-y=1 rather than z-x=2
!!!!!!!!!!!! 1 normalisation
!!!!!!!!!!!! 1 FFT forward
!!!!!!!!!!!! 2 post-processing z-->x direction not pencil
!!!!!!!!!!!! 1 transform z-y
!!!!!!!!!!!! 1 post-processing y
!!!!!!!!!!!! 1 poisson solver
!!!!!!!!!!!! 1 matrice_refinement
!!!!!!!!!!!! 3 inversion5_v1 inversion5_v2
!!!!!!!!!!!! 1 post-processing backward Y
!!!!!!!!!!!! 1 transform y-z
!!!!!!!!!!!! 2 post-processing backward x---->z direction
!!!!!!!!!!!! 1 FFT backward
!!!!!!!!!!!! 1+1 transform |destination: modify the rhs
!!!!!!!!!!!! output real rhs(:,:,:) in the z pencil
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!so the key factor a a2 a3 ax bx
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ay by
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! az bz
call transpose_z_to_y(rhs,rw2,ph)
do k=ph%yst(3),ph%yen(3)
do i=ph%yst(1),ph%yen(1)
do j=1,ny/2
rw2b(i,j,k)=rw2(i,2*(j-1)+1,k)
enddo
do j=ny/2+1,ny
rw2b(i,j,k)=rw2(i,2*ny-2*j+2,k)
enddo
enddo
end do
call transpose_y_to_z(rw2b,rhs,ph)
if (.not. fft_initialised) then
call decomp_2d_fft_init(PHYSICAL_IN_Z,nx,ny,nz)
fft_initialised = .true.
end if
! compute r2c transform
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call decomp_2d_fft_3d(rhs,cw1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! normalisation
cw1 = cw1 / real(nx, kind=mytype) /real(ny, kind=mytype) &
/ real(nz, kind=mytype)
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'START',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! post-processing in spectral space
! POST PROCESSING IN Z
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)+tmp2*az(k), &
tmp2*bz(k)-tmp1*az(k), kind=mytype)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'after z',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! POST PROCESSING IN X
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bx(i)+tmp2*ax(i), &
tmp2*bx(i)-tmp1*ax(i), kind=mytype)
if (i.gt.(nx/2+1)) cw1(i,j,k)=-cw1(i,j,k)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'after x',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! POST PROCESSING IN Y
! NEED TO BE IN Y PENCILS!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!sp!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call transpose_x_to_y(cw1,cw2,sp)
!!!!!!!!!!!!!!!!!!!!!sp!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do k = sp%yst(3), sp%yen(3)
do i = sp%yst(1), sp%yen(1)
cw2b(i,1,k)=cw2(i,1,k)
do j = 2,ny
tmp1 = real(cw2(i,j,k), kind=mytype)
tmp2 = aimag(cw2(i,j,k))
tmp3 = real(cw2(i,ny-j+2,k), kind=mytype)
tmp4 = aimag(cw2(i,ny-j+2,k))
xx1=tmp1*by(j)/2._mytype
xx2=tmp1*ay(j)/2._mytype
xx3=tmp2*by(j)/2._mytype
xx4=tmp2*ay(j)/2._mytype
xx5=tmp3*by(j)/2._mytype
xx6=tmp3*ay(j)/2._mytype
xx7=tmp4*by(j)/2._mytype
xx8=tmp4*ay(j)/2._mytype
cw2b(i,j,k) = cmplx(xx1+xx4+xx5-xx8,-xx2+xx3+xx6+xx7, &
kind=mytype)
end do
end do
end do
#ifdef DEBUG
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
if (abs(cw2b(i,j,k)) > 1.0e-4) then
write(*,100) 'after y',i,j,k,cw2b(i,j,k)
print *,kxyz(i,j,k)
end if
end do
end do
end do
#endif
if (istret==0) then
! Solve Poisson
! doing wave number division in Y-pencil
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
!tmp1=real(zk2(k)+yk2(j)+xk2(i), kind=mytype)
!tmp2=aimag(zk2(k)+yk2(j)+xk2(i))
tmp1=real(kxyz(i,j,k), kind=mytype)
tmp2=aimag(kxyz(i,j,k))
!xyzk=cmplx(tmp1,tmp2, kind=mytype)
!CANNOT DO A DIVISION BY ZERO
if ((abs(tmp1).lt.epsilon).and.(abs(tmp2).lt.epsilon)) then
cw2b(i,j,k)=cmplx(0._mytype,0._mytype, kind=mytype)
end if
if ((abs(tmp1).lt.epsilon).and.(abs(tmp2).ge.epsilon)) then
cw2b(i,j,k)=cmplx(0._mytype, &
aimag(cw2b(i,j,k))/(-tmp2), kind=mytype)
end if
if ((abs(tmp1).ge.epsilon).and.(abs(tmp2).lt.epsilon)) then
cw2b(i,j,k)=cmplx( real(cw2b(i,j,k), kind=mytype) &
/(-tmp1), 0._mytype, kind=mytype)
end if
if ((abs(tmp1).ge.epsilon).and.(abs(tmp2).ge.epsilon)) then
cw2b(i,j,k)=cmplx( real(cw2b(i,j,k), kind=mytype) &
/(-tmp1), &
aimag(cw2b(i,j,k))/(-tmp2), kind=mytype)
end if
end do
end do
end do
else
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call matrice_refinement()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! do k = sp%yst(3), sp%yen(3)
! do j = 1,ny/2
! do i = sp%yst(1), sp%yen(1)
! print *,i,j,k,a(i,j,k,3)
!! if (nrank.le.1) print *,i,j,k,a(i,j,k,3)
!! if (nrank.gt.1) print *,i+4,j,k,a(i,j,k,3)
! enddo
! enddo
! enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!istret !=3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (istret.ne.3) then
cw2(:,:,:)=0.;cw2c(:,:,:)=0.
do k = sp%yst(3), sp%yen(3)
do j = 1,ny/2
do i = sp%yst(1), sp%yen(1)
cw2(i,j,k)=cw2b(i,2*j-1,k)
cw2c(i,j,k)=cw2b(i,2*j,k)
enddo
enddo
enddo
! do k = sp%yst(3), sp%yen(3)
! do j = 1,ny/2
! do i = sp%yst(1), sp%yen(1)
! if (abs(cw2(i,j,k)) > 1.0e-4) then
! write(*,*) 'before IN',i,j,k,cw2(i,j,k)!*2.
!! end if
! end do
! end do
! end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!for the strech grid!!!!!!!!!!!!!!!!!!!!!!!
call inversion5_v1(a,cw2,sp)
call inversion5_v1(a2,cw2c,sp)
! cw2(1,1,1)=cw2(1,1,1)*0.5
! do k = sp%yst(3), sp%yen(3)
! do j = 1,ny/2
! do i = sp%yst(1), sp%yen(1)
! if (abs(cw2c(i,j,k)) > 1.0e-4) then
! write(*,*) 'after IN',i,j,k,cw2c(i,j,k)!*2.
! end if
! end do
! end do
! end do
cw2b(:,:,:)=0.
do k=sp%yst(3), sp%yen(3)
do j=1,ny-1,2
do i=sp%yst(1), sp%yen(1)
cw2b(i,j,k)=cw2(i,(j+1)/2,k)
enddo
enddo
do j=2,ny,2
do i=sp%yst(1), sp%yen(1)
cw2b(i,j,k)=cw2c(i,j/2,k)
enddo
enddo
enddo
!do k=sp%yst(3), sp%yen(3)
!do i=sp%yst(1), sp%yen(1)
! if ((xkx(i)==0).and.(zkz(k)==0)) then
! ! cw2b(i,1,1)=0.
! ! cw2b(i,ny,1)=0.
! endif
!enddo
!enddo
else
do k = sp%yst(3), sp%yen(3)
do j = 1,ny
do i = sp%yst(1), sp%yen(1)
cw2(i,j,k)=cw2b(i,j,k)
enddo
enddo
enddo
call inversion5_v2(a3,cw2,sp)
do k = sp%yst(3), sp%yen(3)
do j = 1,ny
do i = sp%yst(1), sp%yen(1)
cw2b(i,j,k)=cw2(i,j,k)
enddo
enddo
enddo
endif
endif
! print *,nrank, sp%yst(3),sp%yen(3),sp%yst(1),sp%yen(1)
!we are in Y pencil
do k = sp%yst(3), sp%yen(3)
do i = sp%yst(1), sp%yen(1)
if ((i==nx/2+1).and.(k==nz/2+1)) then
cw2b(i,:,k)=0.
endif
enddo
enddo
#ifdef DEBUG
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
if (abs(cw2b(i,j,k)) > 1.0e-4) then
write(*,100) 'AFTER',i,j,k,cw2b(i,j,k)
print *,kxyz(i,j,k)
end if
end do
end do
end do
#endif
! post-processing backward
! POST PROCESSING IN Y
do k = sp%yst(3), sp%yen(3)
do i = sp%yst(1), sp%yen(1)
cw2(i,1,k)=cw2b(i,1,k)
do j = 2,ny
tmp1 = real(cw2b(i,j,k), kind=mytype)
tmp2 = aimag(cw2b(i,j,k))
tmp3 = real(cw2b(i,ny-j+2,k), kind=mytype)
tmp4 = aimag(cw2b(i,ny-j+2,k))
xx1=tmp1*by(j)
xx2=tmp1*ay(j)
xx3=tmp2*by(j)
xx4=tmp2*ay(j)
xx5=tmp3*by(j)
xx6=tmp3*ay(j)
xx7=tmp4*by(j)
xx8=tmp4*ay(j)
cw2(i,j,k) = cmplx(xx1-xx4+xx6+xx7,-(-xx2-xx3+xx5-xx8), &
kind=mytype)
end do
end do
end do
! Back to X-pencil
call transpose_y_to_x(cw2,cw1,sp)
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'AFTER Y',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! POST PROCESSING IN X
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bx(i)-tmp2*ax(i), &
tmp2*bx(i)+tmp1*ax(i), kind=mytype)
if (i.gt.(nx/2+1)) cw1(i,j,k)=-cw1(i,j,k)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'AFTER X',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! POST PROCESSING IN Z
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)-tmp2*az(k), &
tmp2*bz(k)+tmp1*az(k), kind=mytype)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'END',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! compute c2r transform, back to physical space
call decomp_2d_fft_3d(cw1,rhs)
! rhs is in Z-pencil but requires global operations in Y
call transpose_z_to_y(rhs,rw2,ph)
do k=ph%yst(3),ph%yen(3)
do i=ph%yst(1),ph%yen(1)
do j=1,ny/2
rw2b(i,2*j-1,k)=rw2(i,j,k)
enddo
do j=1,ny/2
rw2b(i,2*j,k)=rw2(i,ny-j+1,k)
enddo
enddo
end do
call transpose_y_to_z(rw2b,rhs,ph)
! call decomp_2d_fft_finalize
return
end subroutine poisson_010
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!Neumann!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!11x!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!Neumann Neumann Neumann!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!Neumann Neumann periodic!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solving 3D Poisson equation: Neumann in X, Y; Neumann/periodic in Z
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine poisson_11x(rhs, nclz1)
!!!!!!!!!!!!!!twice real to complex
!!!!!!!!!!! 1: transform rhs-rw2-rw1
!!!!!!!!!!! 2: fft rhs-----rw1
!!!!!!!!!!!!!!twice complex to real
!!!!!!!!!!! 1: transfrom rw1-rw2--rhs
!!!!!!!!!!! 2: fft rw1------rhs
!!!!!!!!!!!!! post-processing operations is under the complex condition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!! input real rhs(:.:,:) in the z pencil
!!!!!!!!!!!! 2+2 transform z--y y->x x-y y-z|destination: modify the rhs
!!!!!!!!!!!! 1 normalisation
!!!!!!!!!!!! 1 FFT forward
!!!!!!!!!!!! 1 post-processing z
!!!!!!!!!!!! 1 transpose x--y
!!!!!!!!!!!! 1 post-processing y
!!!!!!!!!!!! 1 transpose y--x
!!!!!!!!!!!! 1 post-processing x (because he thinks he is in X pencil, Zhaoliang said no)
!!!!!!!!!!!! 1 poisson solver
!!!!!!!!!!!! 1 matrice_refinemento
!!!!!!!!!!!! 1 transpse x--y (for strecting because strecthing is in the y direction)
!!!!!!!!!!!! 1 inversion5_v1
!!!!!!!!!!!! 1 transpose y---x
!!!!!!!!!!!!
!!!!!!!!!!!! 1 post-processing x
!!!!!!!!!!!! 1 transpose x-y
!!!!!!!!!!!! 1 post-processing y
!!!!!!!!!!!! 1 transpose y-x
!!!!!!!!!!!! 1 post-processing z
!!!!!!!!!!!!!
!!!!!!!!!!!! 1 FFT backward
!!!!!!!!!!!! 2+2 transform |destination: modify the rhs
!!!!!!!!!!!! output real rhs(:,:,:) in the z pencil
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
integer, intent(IN) :: nclz1
real(mytype), dimension(:,:,:), intent(INOUT) :: rhs
complex(mytype) :: xyzk
real(mytype) :: tmp1, tmp2, tmp3, tmp4
real(mytype) :: xx1,xx2,xx3,xx4,xx5,xx6,xx7,xx8
integer :: nx,ny,nz, i,j,k
100 format(1x,a8,3I4,2F12.6)
nx = nx_global - 1 ! becasuse (110)
ny = ny_global - 1
if (nclz1==1) then !!!!!the free-slip boundary condition
nz = nz_global - 1
else if (nclz1==0) then
nz = nz_global
end if
if (nclz1==1) then
do j=1,ph%zsz(2)
do i=1,ph%zsz(1)
do k=1,nz/2
rw3(i,j,k)=rhs(i,j,2*(k-1)+1)
end do
do k=nz/2+1,nz
rw3(i,j,k)=rhs(i,j,2*nz-2*k+2)
end do
end do
end do
call transpose_z_to_y(rw3,rw2,ph)
else if (nclz1==0) then
call transpose_z_to_y(rhs,rw2,ph)
end if
do k=ph%yst(3),ph%yen(3)
do i=ph%yst(1),ph%yen(1)
do j=1,ny/2
rw2b(i,j,k)=rw2(i,2*(j-1)+1,k)
end do
do j=ny/2+1,ny
rw2b(i,j,k)=rw2(i,2*ny-2*j+2,k)
end do
end do
end do
! the global operations in X
call transpose_y_to_x(rw2b,rw1,ph)
do k=ph%xst(3),ph%xen(3)
do j=ph%xst(2),ph%xen(2)
do i=1,nx/2
rw1b(i,j,k)=rw1(2*(i-1)+1,j,k)
end do
do i=nx/2+1,nx
rw1b(i,j,k)=rw1(2*nx-2*i+2,j,k)
end do
end do
end do
! back to Z-pencil
call transpose_x_to_y(rw1b,rw2,ph)
call transpose_y_to_z(rw2,rhs,ph)
if (.not. fft_initialised) then
call decomp_2d_fft_init(PHYSICAL_IN_Z,nx,ny,nz)
fft_initialised = .true.
end if
! compute r2c transform
call decomp_2d_fft_3d(rhs,cw1)
! normalisation
cw1 = cw1 / real(nx, kind=mytype) /real(ny, kind=mytype) &
/ real(nz, kind=mytype)
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'START',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! post-processing in spectral space
! POST PROCESSING IN Z
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)+tmp2*az(k), &
tmp2*bz(k)-tmp1*az(k), kind=mytype)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'after z',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!while 1 should be in the corresponding pencils
!!!!!!!!!!!!!!!!!!!!!!!!!!!characteristic 2 : there sholud be tmp1*2*3*4
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! POST PROCESSING IN Y
! WE HAVE TO BE IN Y PENCILS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!No we should be the z pencil!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!there is something wrong in the source code!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!But now we are in the z pencil not in
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! the x pencil!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!so,it should be
!!!!!!!!!!!!!!!!!call transpose_z_to_y rather than transpose_x_to_y
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call transpose_x_to_y(cw1,cw2,sp)
do k = sp%yst(3), sp%yen(3)
do i = sp%yst(1), sp%yen(1)
cw2b(i,1,k)=cw2(i,1,k)
do j = 2,ny
tmp1 = real(cw2(i,j,k), kind=mytype)
tmp2 = aimag(cw2(i,j,k))
tmp3 = real(cw2(i,ny-j+2,k), kind=mytype)
tmp4 = aimag(cw2(i,ny-j+2,k))
xx1=tmp1*by(j)/2._mytype
xx2=tmp1*ay(j)/2._mytype
xx3=tmp2*by(j)/2._mytype
xx4=tmp2*ay(j)/2._mytype
xx5=tmp3*by(j)/2._mytype
xx6=tmp3*ay(j)/2._mytype
xx7=tmp4*by(j)/2._mytype
xx8=tmp4*ay(j)/2._mytype
cw2b(i,j,k) = cmplx(xx1+xx4+xx5-xx8,-xx2+xx3+xx6+xx7, &
kind=mytype)
end do
end do
end do
! back to X-pencil
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!the same wrong ,now you should
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!go back to the z pencil rather than the x pencil
call transpose_y_to_x(cw2b,cw1,sp)
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'after y',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! POST PROCESSING IN X
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
cw1b(1,j,k)=cw1(1,j,k)
do i = 2,nx
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
tmp3 = real(cw1(nx-i+2,j,k), kind=mytype)
tmp4 = aimag(cw1(nx-i+2,j,k))
xx1=tmp1*bx(i)/2._mytype
xx2=tmp1*ax(i)/2._mytype
xx3=tmp2*bx(i)/2._mytype
xx4=tmp2*ax(i)/2._mytype
xx5=tmp3*bx(i)/2._mytype
xx6=tmp3*ax(i)/2._mytype
xx7=tmp4*bx(i)/2._mytype
xx8=tmp4*ax(i)/2._mytype
cw1b(i,j,k) = cmplx(xx1+xx4+xx5-xx8,-xx2+xx3+xx6+xx7, &
kind=mytype)
end do
end do
end do
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1b(i,j,k)) > 1.0e-4) then
write(*,*) 'BEFORE',i,j,k,cw1b(i,j,k)
end if
end do
end do
end do
#endif
if (istret==0) then
! Solve Poisson
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
!tmp1=real(zk2(k)+yk2(j)+xk2(i), kind=mytype)
!tmp2=aimag(zk2(k)+yk2(j)+xk2(i))
tmp1=real(kxyz(i,j,k), kind=mytype)
tmp2=aimag(kxyz(i,j,k))
!xyzk=cmplx(tmp1,tmp2, kind=mytype)
!CANNOT DO A DIVISION BY ZERO
if ((abs(tmp1).lt.epsilon).and.(abs(tmp2).lt.epsilon)) then
cw1b(i,j,k)=cmplx(0._mytype,0._mytype, kind=mytype)
end if
if ((abs(tmp1).lt.epsilon).and.(abs(tmp2).ge.epsilon)) then
cw1b(i,j,k)=cmplx(0._mytype, &
aimag(cw1b(i,j,k))/(-tmp2), kind=mytype)
end if
if ((abs(tmp1).ge.epsilon).and.(abs(tmp2).lt.epsilon)) then
cw1b(i,j,k)=cmplx( real(cw1b(i,j,k), kind=mytype) &
/(-tmp1), 0._mytype, kind=mytype)
end if
if ((abs(tmp1).ge.epsilon).and.(abs(tmp2).ge.epsilon)) then
cw1b(i,j,k)=cmplx( real(cw1b(i,j,k), kind=mytype) &
/(-tmp1), &
aimag(cw1b(i,j,k))/(-tmp2), kind=mytype)
end if
end do
end do
end do
else
call matrice_refinement()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11so we should z->y pencil but not x->ypencil
! the stretching is only working in Y pencils
!!!!!!!!!!!!!!!why the default the pencil is x pencil rather than z pencil!!!
call transpose_x_to_y(cw1b,cw2b,sp)
!we are now in Y pencil
if (istret.ne.3) then
cw2(:,:,:)=0.;cw2c(:,:,:)=0.
do k = sp%yst(3), sp%yen(3)
do j = 1,ny/2
do i = sp%yst(1), sp%yen(1)
cw2(i,j,k)=cw2b(i,2*j-1,k)
cw2c(i,j,k)=cw2b(i,2*j,k)
enddo
enddo
enddo
call inversion5_v1(a,cw2,sp)
call inversion5_v1(a2,cw2c,sp)
cw2b(:,:,:)=0.
do k=sp%yst(3), sp%yen(3)
do j=1,ny-1,2
do i=sp%yst(1), sp%yen(1)
cw2b(i,j,k)=cw2(i,(j+1)/2,k)
enddo
enddo
do j=2,ny,2
do i=sp%yst(1), sp%yen(1)
cw2b(i,j,k)=cw2c(i,j/2,k)
enddo
enddo
enddo
else
cw2(:,:,:)=0.
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
cw2(i,j,k)=cw2b(i,j,k)
enddo
enddo
enddo
call inversion5_v2(a3,cw2,sp)
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
cw2b(i,j,k)=cw2(i,j,k)
enddo
enddo
enddo
endif
!we have to go back in X pencils
call transpose_y_to_x(cw2b,cw1b,sp)
endif
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1b(i,j,k)) > 1.0e-6) then
write(*,*) 'AFTER',i,j,k,cw1b(i,j,k)
end if
end do
end do
end do
#endif
!stop
! post-processing backward
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
cw1(1,j,k)=cw1b(1,j,k)
do i = 2,nx
tmp1 = real(cw1b(i,j,k), kind=mytype)
tmp2 = aimag(cw1b(i,j,k))
tmp3 = real(cw1b(nx-i+2,j,k), kind=mytype)
tmp4 = aimag(cw1b(nx-i+2,j,k))
xx1=tmp1*bx(i)
xx2=tmp1*ax(i)
xx3=tmp2*bx(i)
xx4=tmp2*ax(i)
xx5=tmp3*bx(i)
xx6=tmp3*ax(i)
xx7=tmp4*bx(i)
xx8=tmp4*ax(i)
cw1(i,j,k) = cmplx(xx1-xx4+xx6+xx7,-(-xx2-xx3+xx5-xx8), &
kind=mytype)
end do
end do
end do
#ifdef DEBUG
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
if (abs(cw1(i,j,k)) > 1.0e-4) then
write(*,100) 'AFTER X',i,j,k,cw1(i,j,k)
end if
end do
end do
end do
#endif
! POST PROCESSING IN Y
! NEED to be in Y-pencil
call transpose_x_to_y(cw1,cw2,sp)
do k = sp%yst(3), sp%yen(3)
do i = sp%yst(1), sp%yen(1)
cw2b(i,1,k)=cw2(i,1,k)
do j = 2,ny
tmp1 = real(cw2(i,j,k), kind=mytype)
tmp2 = aimag(cw2(i,j,k))
tmp3 = real(cw2(i,ny-j+2,k), kind=mytype)
tmp4 = aimag(cw2(i,ny-j+2,k))
xx1=tmp1*by(j)
xx2=tmp1*ay(j)
xx3=tmp2*by(j)
xx4=tmp2*ay(j)
xx5=tmp3*by(j)
xx6=tmp3*ay(j)
xx7=tmp4*by(j)
xx8=tmp4*ay(j)
cw2b(i,j,k) = cmplx(xx1-xx4+xx6+xx7,-(-xx2-xx3+xx5-xx8), &
kind=mytype)
end do
end do
end do
#ifdef DEBUG
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
if (abs(cw2b(i,j,k)) > 1.0e-4) then
write(*,100) 'AFTER Y',i,j,k,cw2b(i,j,k)
end if
end do
end do
end do
#endif
! back to X-pencil
call transpose_y_to_x(cw2b,cw1,sp)
! POST PROCESSING IN Z
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
tmp1 = real(cw1(i,j,k), kind=mytype)
tmp2 = aimag(cw1(i,j,k))
cw1(i,j,k) = cmplx(tmp1*bz(k)-tmp2*az(k), &
tmp2*bz(k)+tmp1*az(k), kind=mytype)
#ifdef DEBUG
if (abs(cw1(i,j,k)) > 1.0e-4) &
write(*,100) 'END',i,j,k,cw1(i,j,k)
#endif
end do
end do
end do
! compute c2r transform, back to physical space
call decomp_2d_fft_3d(cw1,rhs)
if (nclz1==1) then
do j=1,ph%zsz(2)
do i=1,ph%zsz(1)
do k=1,nz/2
rw3(i,j,2*k-1)=rhs(i,j,k)
end do
do k=1,nz/2
rw3(i,j,2*k)=rhs(i,j,nz-k+1)
end do
end do
end do
call transpose_z_to_y(rw3,rw2,ph)
else if (nclz1==0) then
call transpose_z_to_y(rhs,rw2,ph)
end if
do k=ph%yst(3),ph%yen(3)
do i=ph%yst(1),ph%yen(1)
do j=1,ny/2
rw2b(i,2*j-1,k)=rw2(i,j,k)
end do
do j=1,ny/2
rw2b(i,2*j,k)=rw2(i,ny-j+1,k)
end do
enddo
end do
call transpose_y_to_x(rw2b,rw1,ph)
do k=ph%xst(3),ph%xen(3)
do j=ph%xst(2),ph%xen(2)
do i=1,nx/2
rw1b(2*i-1,j,k)=rw1(i,j,k)
enddo
do i=1,nx/2
rw1b(2*i,j,k)=rw1(nx-i+1,j,k)
enddo
enddo
end do
call transpose_x_to_y(rw1b,rw2,ph)
call transpose_y_to_z(rw2,rhs,ph)
! call decomp_2d_fft_finalize
return
end subroutine poisson_11x
subroutine abxyz(ax,ay,az,bx,by,bz,nx,ny,nz,bcx,bcy,bcz)
use param
implicit none
integer, intent(IN) :: nx,ny,nz
integer, intent(IN) :: bcx,bcy,bcz
real(mytype), dimension(:), intent(OUT) :: ax,bx
real(mytype), dimension(:), intent(OUT) :: ay,by
real(mytype), dimension(:), intent(OUT) :: az,bz
integer :: i,j,k
if (bcx==0) then
do i=1,nx
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!! generate the x direction coeficiency!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ax(i) = sin(real(i-1, kind=mytype)*PI/real(nx, kind=mytype))
bx(i) = cos(real(i-1, kind=mytype)*PI/real(nx, kind=mytype))
end do
else if (bcx==1) then
do i=1,nx
!!!!!!!!!!!!!!!!!!!!!!!!!one and a half of PI
ax(i) = sin(real(i-1, kind=mytype)*PI/2.0_mytype/ &
real(nx, kind=mytype))
bx(i) = cos(real(i-1, kind=mytype)*PI/2.0_mytype/ &
real(nx, kind=mytype))
end do
end if
if (bcy==0) then
do j=1,ny
ay(j) = sin(real(j-1, kind=mytype)*PI/real(ny, kind=mytype))
by(j) = cos(real(j-1, kind=mytype)*PI/real(ny, kind=mytype))
end do
else if (bcy==1) then
do j=1,ny
ay(j) = sin(real(j-1, kind=mytype)*PI/2.0_mytype/ &
real(ny, kind=mytype))
by(j) = cos(real(j-1, kind=mytype)*PI/2.0_mytype/ &
real(ny, kind=mytype))
end do
end if
if (bcz==0) then
do k=1,nz
az(k) = sin(real(k-1, kind=mytype)*PI/real(nz, kind=mytype))
bz(k) = cos(real(k-1, kind=mytype)*PI/real(nz, kind=mytype))
end do
else if (bcz==1) then
do k=1,nz
az(k) = sin(real(k-1, kind=mytype)*PI/2.0_mytype/ &
real(nz, kind=mytype))
bz(k) = cos(real(k-1, kind=mytype)*PI/2.0_mytype/ &
real(nz, kind=mytype))
end do
end if
return
end subroutine abxyz
! ***********************************************************
!
subroutine waves ()
!
!***********************************************************
USE derivX
USE derivY
USE derivZ
USE param
USE decomp_2d
USE variables
use decomp_2d_fft
implicit none
integer :: i,j,k
real(mytype) :: w,wp,w1,w1p
complex(mytype) :: xyzk
complex(mytype) :: ytt,xtt,ztt,yt1,xt1,yt2,xt2
complex(mytype) :: xtt1,ytt1,ztt1,zt1,zt2,tmp1,tmp2,tmp3
complex(mytype) :: tmp4,tmp5,tmp6
xkx(:)=0. ; xk2(:)=0. ; yky(:)=0. ; yk2(:)=0.
zkz(:)=0. ; zk2(:)=0.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!WAVE NUMBER IN X
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (nclx==0) then
do i=1,nx/2+1
w=2.*pi*(i-1)/nx
wp=acix6*2.*dx*sin(w/2.)+(bcix6*2.*dx)*sin(3./2.*w)
wp=wp/(1.+2.*alcaix6*cos(w))
xkx(i)=cmplx(nx*wp/xlx,nx*wp/xlx, kind=mytype)
exs(i)=cmplx(nx*w/xlx,nx*w/xlx, kind=mytype)
xk2(i)=cmplx((nx*wp/xlx)**2,(nx*wp/xlx)**2, kind=mytype)
enddo
do i=nx/2+2,nx
xkx(i)=xkx(nx-i+2)
exs(i)=exs(nx-i+2)
xk2(i)=xk2(nx-i+2)
enddo
else
do i=1,nx
w=2.*pi*0.5*(i-1)/nxm
wp=acix6*2.*dx*sin(w/2.)+(bcix6*2.*dx)*sin(3./2.*w)
wp=wp/(1.+2.*alcaix6*cos(w))
xkx(i)=cmplx(nxm*wp/xlx,nxm*wp/xlx, kind=mytype)
exs(i)=cmplx(nxm*w/xlx,nxm*w/xlx, kind=mytype)
xk2(i)=cmplx((nxm*wp/xlx)**2,(nxm*wp/xlx)**2, kind=mytype)
enddo
xkx(1)=0.
exs(1)=0.
xk2(1)=0.
endif
!WAVE NUMBER IN Y
if (ncly==0) then
do j=1,ny/2+1
w=2.*pi*(j-1)/ny
wp=aciy6*2.*dy*sin(w/2.)+(bciy6*2.*dy)*sin(3./2.*w)
wp=wp/(1.+2.*alcaiy6*cos(w))
if (istret==0) yky(j)=cmplx(ny*wp/yly,ny*wp/yly, kind=mytype)
if (istret.ne.0) yky(j)=cmplx(ny*wp,ny*wp, kind=mytype)
eys(j)=cmplx(ny*w/yly,ny*w/yly, kind=mytype)
yk2(j)=cmplx((ny*wp/yly)**2,(ny*wp/yly)**2, kind=mytype)
enddo
do j=ny/2+2,ny
yky(j)=yky(ny-j+2)
eys(j)=eys(ny-j+2)
yk2(j)=yk2(ny-j+2)
enddo
else
do j=1,ny
w=2.*pi*0.5*(j-1)/nym
wp=aciy6*2.*dy*sin(w/2.)+(bciy6*2.*dy)*sin(3./2.*w)
wp=wp/(1.+2.*alcaiy6*cos(w))
if (istret==0) yky(j)=cmplx(nym*wp/yly,nym*wp/yly, kind=mytype)
if (istret.ne.0) yky(j)=cmplx(nym*wp,nym*wp, kind=mytype)
eys(j)=cmplx(nym*w/yly,nym*w/yly, kind=mytype)
yk2(j)=cmplx((nym*wp/yly)**2,(nym*wp/yly)**2, kind=mytype)
enddo
yky(1)=0.
eys(1)=0.
yk2(1)=0.
endif
!WAVE NUMBER IN Z
if (nclz==0) then
do k=1,nz/2+1
w=2.*pi*(k-1)/nz
wp=aciz6*2.*dz*sin(w/2.)+(bciz6*2.*dz)*sin(3./2.*w)
wp=wp/(1.+2.*alcaiz6*cos(w))
zkz(k)=cmplx(nz*wp/zlz,nz*wp/zlz, kind=mytype)
ezs(k)=cmplx(nz*w/zlz,nz*w/zlz, kind=mytype)
zk2(k)=cmplx((nz*wp/zlz)**2,(nz*wp/zlz)**2, kind=mytype)
enddo
else
do k=1,nz/2+1
w=2.*pi*0.5*(k-1)/nzm
w1=2.*pi*0.5*(nzm-k+1)/nzm
wp=aciz6*2.*dz*sin(w/2.)+(bciz6*2.*dz)*sin(3./2.*w)
wp=wp/(1.+2.*alcaiz6*cos(w))
w1p=aciz6*2.*dz*sin(w1/2.)+(bciz6*2.*dz)*sin(3./2.*w1)
w1p=w1p/(1.+2.*alcaiz6*cos(w1))
zkz(k)=cmplx(nzm*wp/zlz,-nzm*w1p/zlz, kind=mytype)
ezs(k)=cmplx(nzm*w/zlz,nzm*w1/zlz, kind=mytype)
zk2(k)=cmplx((nzm*wp/zlz)**2,(nzm*w1p/zlz)**2, kind=mytype)
enddo
endif
!
!if (nrank==0) then
! do i=1,nx
! print *,i,ezs(i)
! enddo
!endif
!stop
if ((nclx==0).and.(nclz==0).and.((ncly==1).or.(ncly==2))) then
do k = sp%yst(3), sp%yen(3)
do j = sp%yst(2), sp%yen(2)
do i = sp%yst(1), sp%yen(1)
xtt=cmplx((bicix6*2.*cos(real(exs(i))*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)),&
(bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)), kind=mytype)
ytt=cmplx((biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)),&
(biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)), kind=mytype)
ztt=cmplx((biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)),&
(biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)), kind=mytype)
xtt1=cmplx((aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)),&
(aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)), kind=mytype)
ytt1=cmplx((aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)),&
(aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)), kind=mytype)
ztt1=cmplx((aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)),&
(aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)), kind=mytype)
xt1=cmplx((1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)),&
(1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)), kind=mytype)
yt1=cmplx((1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)),&
(1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)), kind=mytype)
zt1=cmplx((1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)),&
(1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)), kind=mytype)
xt2=xk2(i)*((((ytt1+ytt)/yt1)*((ztt1+ztt)/zt1))**2)
yt2=yk2(j)*((((xtt1+xtt)/xt1)*((ztt1+ztt)/zt1))**2)
zt2=zk2(k)*((((xtt1+xtt)/xt1)*((ytt1+ytt)/yt1))**2)
xyzk=xt2+yt2+zt2
kxyz(i,j,k)=xyzk
! print *,i,j,k, kxyz(i,j,k)
enddo
enddo
enddo
else
if (nclz==0) then
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
xtt=cmplx((bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)),&
(bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)), kind=mytype)
ytt=cmplx((biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)),&
(biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)), kind=mytype)
ztt=cmplx((biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)),&
(biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)), kind=mytype)
xtt1=cmplx((aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)),&
(aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)), kind=mytype)
ytt1=cmplx((aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)),&
(aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)), kind=mytype)
ztt1=cmplx((aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)),&
(aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)), kind=mytype)
xt1=cmplx((1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)),&
(1.+2.*ailcaix6*cos(real(exs(i))*dx)), kind=mytype)
yt1=cmplx((1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)),&
(1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)), kind=mytype)
zt1=cmplx((1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)),&
(1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)), kind=mytype)
xt2=xk2(i)*((((ytt1+ytt)/yt1)*((ztt1+ztt)/zt1))**2)
yt2=yk2(j)*((((xtt1+xtt)/xt1)*((ztt1+ztt)/zt1))**2)
zt2=zk2(k)*((((xtt1+xtt)/xt1)*((ytt1+ytt)/yt1))**2)
xyzk=xt2+yt2+zt2
kxyz(i,j,k)=xyzk
! print *,i,j,k, kxyz(i,j,k)
enddo
enddo
enddo
else
do k = sp%xst(3),sp%xen(3)
do j = sp%xst(2),sp%xen(2)
do i = sp%xst(1),sp%xen(1)
xtt=cmplx((bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)),&
(bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)), kind=mytype)
ytt=cmplx((biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)),&
(biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)), kind=mytype)
!
ztt=cmplx((biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)),&
(biciz6*2.*cos(aimag(ezs(k))*3.*dz/2.)+&
ciciz6*2.*cos(aimag(ezs(k))*5.*dz/2.)), kind=mytype)
!
xtt1=cmplx((aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)),&
(aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)), kind=mytype)
ytt1=cmplx((aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)),&
(aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)), kind=mytype)
!
ztt1=cmplx((aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)),&
(aiciz6*2.*cos(aimag(ezs(k))*dz/2.)), kind=mytype)
!
xt1=cmplx((1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)),&
(1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)), kind=mytype)
yt1=cmplx((1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)),&
(1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)), kind=mytype)
zt1=cmplx((1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)),&
(1.+2.*ailcaiz6*cos(aimag(ezs(k))*dz)), kind=mytype)
tmp1=cmplx(real(ztt1+ztt, kind=mytype)/real(zt1, kind=mytype),&
aimag(ztt1+ztt)/aimag(zt1), kind=mytype)
tmp2=cmplx(real(ytt1+ytt, kind=mytype)/real(yt1, kind=mytype),&
real(ytt1+ytt, kind=mytype)/real(yt1, kind=mytype), kind=mytype)
tmp3=cmplx(real(xtt1+xtt, kind=mytype)/real(xt1, kind=mytype),&
real(xtt1+xtt, kind=mytype)/real(xt1, kind=mytype), kind=mytype)
tmp4=cmplx((real(tmp1, kind=mytype)*real(tmp2, kind=mytype))**2,(aimag(tmp1)*aimag(tmp2))**2, kind=mytype)
tmp5=cmplx((real(tmp1, kind=mytype)*real(tmp3, kind=mytype))**2,(aimag(tmp1)*aimag(tmp3))**2, kind=mytype)
tmp6=cmplx((real(tmp3, kind=mytype)*real(tmp2, kind=mytype))**2,(aimag(tmp3)*aimag(tmp2))**2, kind=mytype)
tmp1=cmplx(real(tmp4, kind=mytype)*real(xk2(i), kind=mytype),aimag(tmp4)*aimag(xk2(i)), kind=mytype)
tmp2=cmplx(real(tmp5, kind=mytype)*real(yk2(j), kind=mytype),aimag(tmp5)*aimag(yk2(j)), kind=mytype)
tmp3=cmplx(real(tmp6, kind=mytype)*real(zk2(k), kind=mytype),aimag(tmp6)*aimag(zk2(k)), kind=mytype)
xyzk=tmp1+tmp2+tmp3
kxyz(i,j,k)=xyzk
! print *,i,j,k,zt1,yt1
enddo
enddo
enddo
endif
endif
! do k=1,1!nz
! do j=1,ny
! do i=1,1!!nx
! print *,j,a(i,j,k,3),kxyz(i,j,k)
! enddo
! enddo
! enddo
end subroutine waves
!**************************************************************************
!
subroutine matrice_refinement()
!
!**************************************************************************
USE decomp_2d
USE variables
USE param
USE var
USE MPI
USE derivX
USE derivY
USE derivZ
implicit none
integer :: i,j,k
complex(mytype),dimension(sp%yst(1):sp%yen(1)) :: transx
complex(mytype),dimension(sp%yst(2):sp%yen(2)) :: transy
complex(mytype),dimension(sp%yst(3):sp%yen(3)) :: transz
real(mytype) :: xa0,xa1
complex(mytype) :: ytt,xtt,ztt,yt1,xt1,yt2,xt2
complex(mytype) :: xtt1,ytt1,ztt1,zt1,zt2,tmp1,tmp2,tmp3
do i = sp%yst(1),sp%yen(1)
xtt=cmplx((bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)),&
(bicix6*2.*cos(real(exs(i), kind=mytype)*3.*dx/2.)+&
cicix6*2.*cos(real(exs(i), kind=mytype)*5.*dx/2.)), kind=mytype)
xtt1=cmplx((aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)),&
(aicix6*2.*cos(real(exs(i), kind=mytype)*dx/2.)), kind=mytype)
xt1=cmplx((1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)),&
(1.+2.*ailcaix6*cos(real(exs(i), kind=mytype)*dx)), kind=mytype)
transx(i)=cmplx(real(xtt1+xtt, kind=mytype)/real(xt1, kind=mytype),&
real(xtt1+xtt, kind=mytype)/real(xt1, kind=mytype), kind=mytype)!(xtt+xtt)/xt1
enddo
do j = sp%yst(2),sp%yen(2)
ytt=cmplx((biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)),&
(biciy6*2.*cos(real(eys(j), kind=mytype)*3.*dy/2.)+&
ciciy6*2.*cos(real(eys(j), kind=mytype)*5.*dy/2.)), kind=mytype)
ytt1=cmplx((aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)),&
(aiciy6*2.*cos(real(eys(j), kind=mytype)*dy/2.)), kind=mytype)
yt1=cmplx((1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)),&
(1.+2.*ailcaiy6*cos(real(eys(j), kind=mytype)*dy)), kind=mytype)
transy(j)=cmplx(real(ytt1+ytt, kind=mytype)/real(yt1, kind=mytype),&
real(ytt1+ytt, kind=mytype)/real(yt1, kind=mytype), kind=mytype)!(ytt+ytt)/yt1
enddo
if (nclz==0) then
do k = sp%yst(3),sp%yen(3)
ztt=cmplx((biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)),&
(biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)), kind=mytype)
ztt1=cmplx((aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)),&
(aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)), kind=mytype)
zt1=cmplx((1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)),&
(1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)), kind=mytype)
transz(k)=cmplx(real(ztt1+ztt, kind=mytype)/real(zt1, kind=mytype),&
aimag(ztt1+ztt)/aimag(zt1), kind=mytype)!(ztt+ztt)/zt1
enddo
else
do k = sp%yst(3),sp%yen(3)
ztt=cmplx((biciz6*2.*cos(real(ezs(k), kind=mytype)*3.*dz/2.)+&
ciciz6*2.*cos(real(ezs(k), kind=mytype)*5.*dz/2.)),&
(biciz6*2.*cos(aimag(ezs(k))*3.*dz/2.)+&
ciciz6*2.*cos(aimag(ezs(k))*5.*dz/2.)), kind=mytype)
ztt1=cmplx((aiciz6*2.*cos(real(ezs(k), kind=mytype)*dz/2.)),&
(aiciz6*2.*cos(aimag(ezs(k))*dz/2.)), kind=mytype)
zt1=cmplx((1.+2.*ailcaiz6*cos(real(ezs(k), kind=mytype)*dz)),&
(1.+2.*ailcaiz6*cos(aimag(ezs(k))*dz)), kind=mytype)
transz(k)=cmplx(real(ztt1+ztt, kind=mytype)/real(zt1, kind=mytype),&
aimag(ztt1+ztt)/aimag(zt1), kind=mytype)!(ztt+ztt)/zt1
enddo
endif
if ((istret==1).or.(istret==2)) then
xa0=alpha/pi+1./2./beta/pi
if (istret==1) xa1=1./4./beta/pi
if (istret==2) xa1=-1./4./beta/pi
!
!construction of the pentadiagonal matrice
!
do k=sp%yst(3),sp%yen(3)
do j=1,ny/2
do i=sp%yst(1),sp%yen(1)
cw22(i,j,k)=cmplx(real(yky(2*j-1), kind=mytype)*real(transx(i), kind=mytype)*real(transz(k), kind=mytype),&
aimag(yky(2*j-1))*aimag(transx(i))*aimag(transz(k)), kind=mytype)
cw2(i,j,k)=cmplx(real(yky(2*j), kind=mytype)*real(transx(i), kind=mytype)*real(transz(k), kind=mytype),&
aimag(yky(2*j))*aimag(transx(i))*aimag(transz(k)), kind=mytype)
enddo
enddo
enddo
!main diagonal
do k=sp%yst(3),sp%yen(3)
do j=2,ny/2-1
do i=sp%yst(1),sp%yen(1)
a(i,j,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(2*j-1), kind=mytype)*real(transy(2*j-1), kind=mytype)&
*real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(2*j-1), kind=mytype)*real(transy(2*j-1), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw22(i,j,k), kind=mytype)*real(cw22(i,j,k), kind=mytype)*xa0*xa0-&
xa1*xa1*(real(cw22(i,j,k), kind=mytype)*real(cw22(i,j-1,k), kind=mytype)+real(cw22(i,j,k), kind=mytype)*&
real(cw22(i,j+1,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(2*j-1))*aimag(transy(2*j-1))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(2*j-1))*aimag(transy(2*j-1))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw22(i,j,k))*aimag(cw22(i,j,k))*xa0*xa0-&
xa1*xa1*(aimag(cw22(i,j,k))*aimag(cw22(i,j-1,k))+aimag(cw22(i,j,k))*aimag(cw22(i,j+1,k))), kind=mytype)
a2(i,j,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(2*j), kind=mytype)*real(transy(2*j), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(2*j), kind=mytype)*real(transy(2*j), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw2(i,j,k), kind=mytype)*real(cw2(i,j,k), kind=mytype)*xa0*xa0-&
xa1*xa1*(real(cw2(i,j,k), kind=mytype)*real(cw2(i,j-1,k), kind=mytype)+real(cw2(i,j,k), kind=mytype)*&
real(cw2(i,j+1,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(2*j))*aimag(transy(2*j))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(2*j))*aimag(transy(2*j))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw2(i,j,k))*aimag(cw2(i,j,k))*xa0*xa0-&
xa1*xa1*(aimag(cw2(i,j,k))*aimag(cw2(i,j-1,k))+aimag(cw2(i,j,k))*aimag(cw2(i,j+1,k))), kind=mytype)
enddo
enddo
do i=sp%yst(1),sp%yen(1)
a(i,1,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(1), kind=mytype)*real(transy(1), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(1), kind=mytype)*real(transy(1), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw22(i,1,k), kind=mytype)*real(cw22(i,1,k), kind=mytype)*xa0*xa0-xa1*xa1*(real(cw22(i,1,k), kind=mytype)*&
real(cw22(i,2,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(1))*aimag(transy(1))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(1))*aimag(transy(1))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw22(i,1,k))*aimag(cw22(i,1,k))*xa0*xa0-xa1*xa1*(aimag(cw22(i,1,k))*aimag(cw22(i,2,k))), kind=mytype)
a(i,ny/2,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(ny-2), kind=mytype)*real(transy(ny-2), kind=mytype)&
*real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(ny-2), kind=mytype)*real(transy(ny-2), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw22(i,ny/2,k), kind=mytype)*real(cw22(i,ny/2,k), kind=mytype)*xa0*xa0-&
xa1*xa1*(real(cw22(i,ny/2,k), kind=mytype)*real(cw22(i,ny/2-1,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(ny-2))*aimag(transy(ny-2))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(ny-2))*aimag(transy(ny-2))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw22(i,ny/2,k))*aimag(cw22(i,ny/2,k))*xa0*xa0-&
xa1*xa1*(aimag(cw22(i,ny/2,k))*aimag(cw22(i,ny/2-1,k))), kind=mytype)
a2(i,1,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(2), kind=mytype)*real(transy(2), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(2), kind=mytype)*real(transy(2), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw2(i,1,k), kind=mytype)*real(cw2(i,1,k), kind=mytype)*(xa0-xa1)*(xa0+xa1)-xa1*xa1*(real(cw2(i,1,k), kind=mytype)*&
real(cw2(i,2,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(2))*aimag(transy(2))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(2))*aimag(transy(2))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw2(i,1,k))*aimag(cw2(i,1,k))*(xa0-xa1)*(xa0+xa1)-xa1*xa1*(aimag(cw2(i,1,k))*aimag(cw2(i,2,k))), kind=mytype)
a2(i,ny/2,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(ny-1), kind=mytype)*real(transy(ny-1), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(ny-1), kind=mytype)*real(transy(ny-1), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw2(i,ny/2,k), kind=mytype)*real(cw2(i,ny/2,k), kind=mytype)*(xa0+xa1)*(xa0+xa1)-&
xa1*xa1*(real(cw2(i,ny/2,k), kind=mytype)*real(cw2(i,ny/2-1,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(ny-1))*aimag(transy(ny-1))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(ny-1))*aimag(transy(ny-1))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw2(i,ny/2,k))*aimag(cw2(i,ny/2,k))*(xa0+xa1)*(xa0+xa1)-&
xa1*xa1*(aimag(cw2(i,ny/2,k))*aimag(cw2(i,ny/2-1,k))), kind=mytype)
enddo
enddo
!sup diag +1
do k=sp%yst(3),sp%yen(3)
do j=2,ny/2-1
do i=sp%yst(1),sp%yen(1)
a(i,j,k,4)=cmplx(xa0*xa1*(real(cw22(i,j,k), kind=mytype)*real(cw22(i,j+1,k), kind=mytype)+real(cw22(i,j+1,k), kind=mytype)*&
real(cw22(i,j+1,k), kind=mytype)),&
xa0*xa1*(aimag(cw22(i,j,k))*aimag(cw22(i,j+1,k))+aimag(cw22(i,j+1,k))*aimag(cw22(i,j+1,k))), kind=mytype)
a2(i,j,k,4)=cmplx(xa0*xa1*(real(cw2(i,j,k), kind=mytype)*real(cw2(i,j+1,k), kind=mytype)+real(cw2(i,j+1,k), kind=mytype)*&
real(cw2(i,j+1,k), kind=mytype)),&
xa0*xa1*(aimag(cw2(i,j,k))*aimag(cw2(i,j+1,k))+aimag(cw2(i,j+1,k))*aimag(cw2(i,j+1,k))), kind=mytype)
enddo
enddo
do i=sp%yst(1),sp%yen(1)
a(i,1,k,4)=2.*cmplx((xa0*xa1*(real(cw22(i,1,k), kind=mytype)*real(cw22(i,2,k), kind=mytype)+real(cw22(i,2,k), kind=mytype)*&
real(cw22(i,2,k), kind=mytype))),&
(xa0*xa1*(aimag(cw22(i,1,k))*aimag(cw22(i,2,k))+aimag(cw22(i,2,k))*aimag(cw22(i,2,k)))), kind=mytype)
a2(i,1,k,4)=cmplx((xa0-xa1)*xa1*(real(cw2(i,1,k), kind=mytype)*real(cw2(i,2,k), kind=mytype))&
+xa0*xa1*(real(cw2(i,2,k), kind=mytype)*real(cw2(i,2,k), kind=mytype)),&
(xa0-xa1)*xa1*(aimag(cw2(i,1,k))*aimag(cw2(i,2,k)))&
+xa0*xa1*(aimag(cw2(i,2,k))*aimag(cw2(i,2,k))), kind=mytype)
a2(i,ny/2-1,k,4)=cmplx(xa0*xa1*(real(cw2(i,ny/2-1,k), kind=mytype)*real(cw2(i,ny/2,k), kind=mytype))&
+(xa0+xa1)*xa1*(real(cw2(i,ny/2,k), kind=mytype)*real(cw2(i,ny/2,k), kind=mytype)),&
xa0*xa1*(aimag(cw2(i,ny/2-1,k))*aimag(cw2(i,ny/2,k)))&
+(xa0+xa1)*xa1*(aimag(cw2(i,ny/2,k))*aimag(cw2(i,ny/2,k))), kind=mytype)
a2(i,ny/2,k,4)=0.
enddo
enddo
!sup diag +2
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=1,ny/2-2
a(i,j,k,5)=cmplx(-real(cw22(i,j+1,k), kind=mytype)*real(cw22(i,j+2,k), kind=mytype)*xa1*xa1,&
-aimag(cw22(i,j+1,k))*aimag(cw22(i,j+2,k))*xa1*xa1, kind=mytype)
a2(i,j,k,5)=cmplx(-real(cw2(i,j+1,k), kind=mytype)*real(cw2(i,j+2,k), kind=mytype)*xa1*xa1,&
-aimag(cw2(i,j+1,k))*aimag(cw2(i,j+2,k))*xa1*xa1, kind=mytype)
enddo
a(i,1,k,5)=cmplx(real(a(i,1,k,5), kind=mytype)*2.,aimag(a(i,1,k,5))*2., kind=mytype)
a(i,ny/2-1,k,5)=0.
a(i,ny/2,k,5)=0.
a2(i,ny/2-1,k,5)=0.
a2(i,ny/2,k,5)=0.
enddo
enddo
!inf diag -1
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=2,ny/2
a(i,j,k,2)=cmplx(xa0*xa1*(real(cw22(i,j,k), kind=mytype)*real(cw22(i,j-1,k), kind=mytype)+real(cw22(i,j-1,k), kind=mytype)*&
real(cw22(i,j-1,k), kind=mytype)),&
xa0*xa1*(aimag(cw22(i,j,k))*aimag(cw22(i,j-1,k))+aimag(cw22(i,j-1,k))*aimag(cw22(i,j-1,k))), kind=mytype)
a2(i,j,k,2)=cmplx(xa0*xa1*(real(cw2(i,j,k), kind=mytype)*real(cw2(i,j-1,k), kind=mytype)+real(cw2(i,j-1,k), kind=mytype)*&
real(cw2(i,j-1,k), kind=mytype)),&
xa0*xa1*(aimag(cw2(i,j,k))*aimag(cw2(i,j-1,k))+aimag(cw2(i,j-1,k))*aimag(cw2(i,j-1,k))), kind=mytype)
enddo
a(i,1,k,2)=0.
a2(i,1,k,2)=0.
a2(i,2,k,2)=cmplx(xa0*xa1*(real(cw2(i,2,k), kind=mytype)*real(cw2(i,1,k), kind=mytype))&
+(xa0+xa1)*xa1*(real(cw2(i,1,k), kind=mytype)*real(cw2(i,1,k), kind=mytype)),&
xa0*xa1*(aimag(cw2(i,2,k))*aimag(cw2(i,1,k)))&
+(xa0+xa1)*xa1*(aimag(cw2(i,1,k))*aimag(cw2(i,1,k))), kind=mytype)
a2(i,ny/2,k,2)=cmplx((xa0+xa1)*xa1*(real(cw2(i,ny/2,k), kind=mytype)*real(cw2(i,ny/2-1,k), kind=mytype))&
+xa0*xa1*(real(cw2(i,ny/2-1,k), kind=mytype)*real(cw2(i,ny/2-1,k), kind=mytype)),&
(xa0+xa1)*xa1*(aimag(cw2(i,ny/2,k))*aimag(cw2(i,ny/2-1,k)))&
+xa0*xa1*(aimag(cw2(i,ny/2-1,k))*aimag(cw2(i,ny/2-1,k))), kind=mytype)
enddo
enddo
!inf diag -2
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=3,ny/2
a(i,j,k,1)=cmplx(-real(cw22(i,j-1,k), kind=mytype)*real(cw22(i,j-2,k), kind=mytype)*xa1*xa1,&
-aimag(cw22(i,j-1,k))*aimag(cw22(i,j-2,k))*xa1*xa1, kind=mytype)
a2(i,j,k,1)=cmplx(-real(cw2(i,j-1,k), kind=mytype)*real(cw2(i,j-2,k), kind=mytype)*xa1*xa1,&
-aimag(cw2(i,j-1,k))*aimag(cw2(i,j-2,k))*xa1*xa1, kind=mytype)
enddo
a(i,1,k,1)=0.
a(i,2,k,1)=0.
a2(i,1,k,1)=0.
a2(i,2,k,1)=0.
enddo
enddo
!not to have a singular matrice
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
if ((real(xk2(i), kind=mytype)==0.).and.(real(zk2(k), kind=mytype)==0)) then
a(i,1,k,3)=cmplx(1.,1., kind=mytype)
a(i,1,k,4)=0.
a(i,1,k,5)=0.
endif
enddo
enddo
else
xa0=alpha/pi+1./2./beta/pi
xa1=-1./4./beta/pi
!
!construction of the pentadiagonal matrice
!
do k=sp%yst(3),sp%yen(3)
do j=1,nym
do i=sp%yst(1),sp%yen(1)
cw22(i,j,k)=cmplx(real(yky(j), kind=mytype)*real(transx(i), kind=mytype)*real(transz(k), kind=mytype),&
aimag(yky(j))*aimag(transx(i))*aimag(transz(k)), kind=mytype)
enddo
enddo
enddo
!main diagonal
do k=sp%yst(3),sp%yen(3)
do j=2,nym-1
do i=sp%yst(1),sp%yen(1)
a3(i,j,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(j), kind=mytype)*real(transy(j), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(j), kind=mytype)*real(transy(j), kind=mytype)*&
real(transx(i), kind=mytype)*real(transx(i), kind=mytype))&
-real(cw22(i,j,k), kind=mytype)*real(cw22(i,j,k), kind=mytype)*xa0*xa0-&
xa1*xa1*(real(cw22(i,j,k), kind=mytype)*real(cw22(i,j-1,k), kind=mytype)+real(cw22(i,j,k), kind=mytype)*&
real(cw22(i,j+1,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(j))*aimag(transy(j))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(j))*aimag(transy(j))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw22(i,j,k))*aimag(cw22(i,j,k))*xa0*xa0-&
xa1*xa1*(aimag(cw22(i,j,k))*aimag(cw22(i,j-1,k))+aimag(cw22(i,j,k))*aimag(cw22(i,j+1,k))), kind=mytype)
enddo
enddo
enddo
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
a3(i,1,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(1), kind=mytype)*real(transy(1), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(1), kind=mytype)*real(transy(1), kind=mytype)*real(transx(i), kind=mytype)*&
real(transx(i), kind=mytype))&
-real(cw22(i,1,k), kind=mytype)*real(cw22(i,1,k), kind=mytype)*xa0*xa0-xa1*xa1*(real(cw22(i,1,k), kind=mytype)*&
real(cw22(i,2,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(1))*aimag(transy(1))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(1))*aimag(transy(1))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw22(i,1,k))*aimag(cw22(i,1,k))*xa0*xa0-xa1*xa1*(aimag(cw22(i,1,k))*aimag(cw22(i,2,k))), kind=mytype)
a3(i,nym,k,3)=cmplx(-(real(xk2(i), kind=mytype)*real(transy(nym), kind=mytype)*real(transy(nym), kind=mytype)*&
real(transz(k), kind=mytype)*real(transz(k), kind=mytype))&
-(real(zk2(k), kind=mytype)*real(transy(nym), kind=mytype)*real(transy(nym), kind=mytype)*real(transx(i), kind=mytype)*&
real(transx(i), kind=mytype))&
-real(cw22(i,nym,k), kind=mytype)*real(cw22(i,nym,k), kind=mytype)*xa0*xa0-&
xa1*xa1*(real(cw22(i,nym,k), kind=mytype)*real(cw22(i,nym-1,k), kind=mytype)),&
-(aimag(xk2(i))*aimag(transy(nym))*aimag(transy(nym))*aimag(transz(k))*aimag(transz(k)))&
-(aimag(zk2(k))*aimag(transy(nym))*aimag(transy(nym))*aimag(transx(i))*aimag(transx(i)))&
-aimag(cw22(i,nym,k))*aimag(cw22(i,nym,k))*xa0*xa0-&
xa1*xa1*(aimag(cw22(i,nym,k))*aimag(cw22(i,nym-1,k))), kind=mytype)
enddo
enddo
!sup diag +1
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=2,nym-1
a3(i,j,k,4)=cmplx(xa0*xa1*(real(cw22(i,j,k), kind=mytype)*real(cw22(i,j+1,k), kind=mytype)+real(cw22(i,j+1,k), kind=mytype)*&
real(cw22(i,j+1,k), kind=mytype)),&
xa0*xa1*(aimag(cw22(i,j,k))*aimag(cw22(i,j+1,k))+aimag(cw22(i,j+1,k))*aimag(cw22(i,j+1,k))), kind=mytype)
enddo
a3(i,1,k,4)=cmplx((xa0*xa1*(real(cw22(i,1,k), kind=mytype)*real(cw22(i,2,k), kind=mytype)+real(cw22(i,2,k), kind=mytype)*&
real(cw22(i,2,k), kind=mytype))),&
(xa0*xa1*(aimag(cw22(i,1,k))*aimag(cw22(i,2,k))+aimag(cw22(i,2,k))*aimag(cw22(i,2,k)))), kind=mytype)
enddo
enddo
!sup diag +2
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=1,nym-2
a3(i,j,k,5)=cmplx(-real(cw22(i,j+1,k), kind=mytype)*real(cw22(i,j+2,k), kind=mytype)*xa1*xa1,&
-aimag(cw22(i,j+1,k))*aimag(cw22(i,j+2,k))*xa1*xa1, kind=mytype)
enddo
!a3(i,1,k,5)=a3(i,1,k,5)*2.
!a3(i,1,k,5)=0.
a3(i,nym-1,k,5)=0.
a3(i,nym,k,5)=0.
enddo
enddo
!inf diag -1
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=2,nym
a3(i,j,k,2)=cmplx(xa0*xa1*(real(cw22(i,j,k), kind=mytype)*real(cw22(i,j-1,k), kind=mytype)+real(cw22(i,j-1,k), kind=mytype)*&
real(cw22(i,j-1,k), kind=mytype)),&
xa0*xa1*(aimag(cw22(i,j,k))*aimag(cw22(i,j-1,k))+aimag(cw22(i,j-1,k))*aimag(cw22(i,j-1,k))), kind=mytype)
enddo
a3(i,1,k,2)=0.
enddo
enddo
!inf diag -2
do k=sp%yst(3),sp%yen(3)
do i=sp%yst(1),sp%yen(1)
do j=3,nym
a3(i,j,k,1)=cmplx(-real(cw22(i,j-1,k), kind=mytype)*real(cw22(i,j-2,k), kind=mytype)*xa1*xa1,&
-aimag(cw22(i,j-1,k))*aimag(cw22(i,j-2,k))*xa1*xa1, kind=mytype)
enddo
a3(i,1,k,1)=0.
a3(i,2,k,1)=0.
enddo
enddo
!not to have a singular matrice
! do k=sp%yst(3),sp%yen(3)
! do i=sp%yst(1),sp%yen(1)
! if ((xkx(i)==0.).and.(zkz(k)==0)) then
if (nrank==0) then
a3(1,1,1,3)=cmplx(1.,1., kind=mytype)
a3(1,1,1,4)=0.
a3(1,1,1,5)=0.
endif
! endif
! enddo
! enddo
endif
return
end subroutine matrice_refinement
end module decomp_2d_poisson
令狐冲
Engineer of offshore wind turbine technique research
My research interests include distributed energy, wind turbine power generation technique , Computational fluid dynamic and programmable matter.