Temperton circulant tridiagonal algorithm
This commit is contained in:
ignis
parent
56fc50e9f0
commit
3ee83a9fff
2 changed files with 429 additions and 3 deletions
428
m_compact.f90
428
m_compact.f90
|
|
@ -8,10 +8,13 @@
|
|||
lyf,lys,wyf,wys, &
|
||||
lzf,lzs,wzf,wzs
|
||||
! lyzf,lyzs,wyzf,wyzs
|
||||
|
||||
REAL*8, DIMENSION(:), ALLOCATABLE :: zz1, zz2, vz1, vz2
|
||||
|
||||
INTEGER :: nxc,nyc,nzc
|
||||
REAL*8, PARAMETER :: ezero = 1.0e-14
|
||||
|
||||
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp
|
||||
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp,pdfp,pd2fp
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
@ -68,6 +71,22 @@
|
|||
call nonp_lud(3,nzc)
|
||||
ENDIF
|
||||
|
||||
ALLOCATE(vz1(nzc),STAT=ierr)
|
||||
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
||||
|
||||
ALLOCATE(zz1(nzc),STAT=ierr)
|
||||
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
||||
|
||||
ALLOCATE(vz2(nzc),STAT=ierr)
|
||||
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
||||
|
||||
ALLOCATE(zz2(nzc),STAT=ierr)
|
||||
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
|
||||
|
||||
call ctdlu (3.0, nzc, vz1, zz1)
|
||||
|
||||
call ctdlu (11.0/2.0, nzc, vz2, zz2)
|
||||
|
||||
! CALL x_lud
|
||||
! CALL yz_lud
|
||||
|
||||
|
|
@ -144,6 +163,232 @@
|
|||
l(n)=1./d
|
||||
END SUBROUTINE ptdlu
|
||||
|
||||
SUBROUTINE ctdlu(lambda,n,v,z)
|
||||
INTEGER, INTENT(IN) :: n
|
||||
REAL*8, INTENT(IN) :: lambda
|
||||
REAL*8, INTENT(OUT) :: v(n)
|
||||
REAL*8, INTENT(OUT) :: z(n)
|
||||
INTEGER :: i
|
||||
REAL*8 :: sigma, alpha
|
||||
|
||||
v(1)=0.0
|
||||
v(2)=lambda
|
||||
DO i=3,n
|
||||
v(i)=lambda - 1./v(i-1)
|
||||
ENDDO
|
||||
|
||||
if (lambda > 0.0) then
|
||||
alpha = (-lambda + sqrt(lambda**2 - 4.0))/2.0
|
||||
else
|
||||
alpha = (-lambda - sqrt(lambda**2 - 4.0))/2.0
|
||||
endif
|
||||
|
||||
sigma = (1.0 + alpha**2) / (lambda * (1.0 - alpha**2) * (1.0 - alpha**n))
|
||||
|
||||
DO i=1,n
|
||||
z(i) = sigma * (alpha**(i-1) + alpha**(n-i+1))
|
||||
ENDDO
|
||||
END SUBROUTINE ctdlu
|
||||
|
||||
|
||||
SUBROUTINE pdfp(x,xu,xl,dx,h,nd,n,nall,dir)
|
||||
|
||||
INTEGER,INTENT(IN) :: nall,n,nd,dir
|
||||
REAL*8,INTENT(IN) :: h
|
||||
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
|
||||
REAL*8,INTENT(OUT),DIMENSION(nd,2) :: xu, xl
|
||||
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
|
||||
|
||||
INTEGER :: i,j
|
||||
REAL*8 :: r1,r2,h1
|
||||
|
||||
h1=1./h
|
||||
r1=7./3.
|
||||
r2=1./12.
|
||||
|
||||
!if (myid.eq.0) write(*,*) "parallel dfp commuication"
|
||||
if (myid.eq.master) then
|
||||
|
||||
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
elseif (myid.eq.numprocs-1) then
|
||||
|
||||
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xu, 2*nd, MPI_REAL8, 0, 0, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
else
|
||||
|
||||
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
endif
|
||||
|
||||
DO j=1,nd
|
||||
dx(j,1) =(r1*( x(j,2)-xl(j,2)) +r2*( x(j,3)-xl(j,1))) *h1
|
||||
ENDDO
|
||||
|
||||
DO j=1,nd
|
||||
dx(j,2) =(r1*( x(j,3)- x(j,1)) +r2*( x(j,4)-xl(j,2))) *h1
|
||||
ENDDO
|
||||
|
||||
DO j=1,nd
|
||||
dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
|
||||
ENDDO
|
||||
|
||||
DO j=1,nd
|
||||
dx(j,n) =(r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
|
||||
ENDDO
|
||||
|
||||
DO i=3,n-2
|
||||
DO j=1,nd
|
||||
dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))*h1
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
|
||||
IF (dir.eq.1) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
|
||||
IF (dir.eq.2) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
|
||||
IF (dir.eq.3) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
|
||||
|
||||
END SUBROUTINE pdfp
|
||||
|
||||
|
||||
SUBROUTINE pd2fp(x,xu,xl,dx,h,nd,n,nall,dir)
|
||||
INTEGER,INTENT(IN) :: nall,n,nd,dir
|
||||
REAL*8,INTENT(IN) :: h
|
||||
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
|
||||
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
|
||||
INTEGER :: i,j
|
||||
REAL*8 :: h2,r1,r2,t1,t2
|
||||
|
||||
REAL*8,DIMENSION(nd,2) :: xu, xl
|
||||
|
||||
h2=1./(h*h)
|
||||
r1=6.
|
||||
r2=3./8.
|
||||
|
||||
!if (myid.eq.0) write(*,*) "parallel d2fp commuication"
|
||||
if (myid.eq.master) then
|
||||
|
||||
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
elseif (myid.eq.numprocs-1) then
|
||||
|
||||
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xu, 2*nd, MPI_REAL8, 0, 0, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
else
|
||||
|
||||
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
|
||||
MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
|
||||
MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
endif
|
||||
|
||||
|
||||
DO j=1,nd
|
||||
|
||||
t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
|
||||
t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
|
||||
dx(j,1)=(r1*t1+r2*t2)*h2
|
||||
|
||||
ENDDO
|
||||
|
||||
DO j=1,nd
|
||||
|
||||
t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
|
||||
t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
|
||||
dx(j,2)=(r1*t1+r2*t2)*h2
|
||||
|
||||
ENDDO
|
||||
|
||||
DO j=1,nd
|
||||
|
||||
t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2))
|
||||
t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
|
||||
dx(j,n-1)=(r1*t1+r2*t2)*h2
|
||||
|
||||
ENDDO
|
||||
|
||||
DO j=1,nd
|
||||
|
||||
t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
|
||||
t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
|
||||
dx(j,n)=(r1*t1+r2*t2)*h2
|
||||
|
||||
ENDDO
|
||||
|
||||
DO i=3,n-2
|
||||
DO j=1,nd
|
||||
|
||||
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
|
||||
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
|
||||
dx(j,i)=(r1*t1+r2*t2)*h2
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
IF (dir.eq.1) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! x-direction
|
||||
IF (dir.eq.2) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! y-direction
|
||||
IF (dir.eq.3) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! z-direction
|
||||
END SUBROUTINE pd2fp
|
||||
|
||||
|
||||
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
|
||||
INTEGER,INTENT(IN) :: n,nd,dir
|
||||
REAL*8,INTENT(IN) :: h
|
||||
|
|
@ -329,6 +574,187 @@
|
|||
ENDDO
|
||||
END SUBROUTINE ptdslv
|
||||
|
||||
SUBROUTINE ctdslv(r,va,za,nd,n,nall)
|
||||
INTEGER,PARAMETER :: nb = 256
|
||||
INTEGER,INTENT(IN) :: n,nd,nall
|
||||
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
REAL*8,INTENT(IN),DIMENSION(nall) :: va
|
||||
REAL*8,INTENT(IN),DIMENSION(nall) :: za
|
||||
REAL*8,DIMENSION(n) :: v,z
|
||||
INTEGER i,j
|
||||
INTEGER ii,jj
|
||||
REAL*8, DIMENSION(nd) :: sum
|
||||
|
||||
REAL*8, DIMENSION(nd) :: x1, x1_
|
||||
|
||||
INTEGER npart, nbase, nlow, nupp
|
||||
INTEGER :: pid, np
|
||||
REAL*8, DIMENSION(nb) :: r0, r1, sum0
|
||||
REAL*8, DIMENSION(nb,2) :: buf
|
||||
|
||||
pid = myid
|
||||
np = numprocs
|
||||
|
||||
npart = nall / np
|
||||
nbase = pid * npart
|
||||
nlow = pid * npart + 1
|
||||
nupp = (pid + 1) * npart
|
||||
|
||||
v = va(nlow:nupp)
|
||||
z = za(nlow:nupp)
|
||||
|
||||
if (npart.lt.4) then
|
||||
! assertion fail
|
||||
endif
|
||||
|
||||
x1 = 0.0
|
||||
x1_ = 0.0
|
||||
|
||||
do i=1,n
|
||||
x1_ = x1_ + z(i)*r(:,i)
|
||||
enddo
|
||||
|
||||
! call MPI_ALLREDUCE(x1_, x1, nd, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
|
||||
call MPI_REDUCE(x1_, x1, nd, MPI_REAL8, MPI_SUM, master, MPI_COMM_TASK, mpi_err)
|
||||
|
||||
! first process
|
||||
if (myid.eq.0) then
|
||||
|
||||
CALL MPI_ISEND(x1, nd, MPI_REAL8, np-1, 100, MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
|
||||
r(:,1) = x1
|
||||
r(:,2) = r(:,2) - x1
|
||||
|
||||
! last process
|
||||
elseif (pid.eq.(np-1)) then
|
||||
|
||||
CALL MPI_RECV(x1, nd, MPI_REAL8, master, 100, MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
r(:,n) = r(:,n) - x1
|
||||
|
||||
! intermediate process
|
||||
else
|
||||
|
||||
endif
|
||||
|
||||
! first process
|
||||
if (myid.eq.0) then
|
||||
|
||||
DO jj=1,nd,nb
|
||||
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,2)=r(j,2)/v(2)
|
||||
ENDDO
|
||||
|
||||
DO i=3,n
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=(r(j,i)-r(j,i-1))/v(i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
r0 = r(jj:jj+nb-1,n)
|
||||
CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)
|
||||
|
||||
ENDDO
|
||||
|
||||
DO jj=1,nd,nb
|
||||
|
||||
CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
i=n
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i) = r(j,i) - r1(j-jj+1)/v(i)
|
||||
ENDDO
|
||||
|
||||
DO i=n-1,2,-1
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i) = r(j,i) - r(j,i+1)/v(i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
|
||||
! last process
|
||||
elseif (pid.eq.(np-1)) then
|
||||
|
||||
DO jj=1,nd,nb
|
||||
|
||||
CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
i=1
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=(r(j,i) - r1(j-jj+1))/v(i)
|
||||
ENDDO
|
||||
|
||||
DO i=2,n
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=(r(j,i) - r(j,i-1))/v(i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
|
||||
DO jj=1,nd,nb
|
||||
|
||||
DO i=n-1,1,-1
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=r(j,i)-r(j,i+1)/v(i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
r0 = r(jj:jj+nb-1,1)
|
||||
CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)
|
||||
|
||||
ENDDO
|
||||
|
||||
! intermediate process
|
||||
else
|
||||
|
||||
DO jj=1,nd,nb
|
||||
|
||||
CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
i=1
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=(r(j,i)-r1(j-jj+1))/v(i)
|
||||
ENDDO
|
||||
|
||||
DO i=2,n
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=(r(j,i)-r(j,i-1))/v(i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
r0 = r(jj:jj+nb-1,n)
|
||||
CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)
|
||||
|
||||
ENDDO
|
||||
|
||||
DO jj=1,nd,nb
|
||||
|
||||
CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
|
||||
i=n
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=r(j,i)-r1(j-jj+1)/v(i)
|
||||
ENDDO
|
||||
|
||||
DO i=n-1,1,-1
|
||||
DO j=jj,jj+nb-1
|
||||
r(j,i)=r(j,i)-r(j,i+1)/v(i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
r0 = r(jj:jj+nb-1,1)
|
||||
CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)
|
||||
|
||||
ENDDO
|
||||
|
||||
endif
|
||||
|
||||
END SUBROUTINE ctdslv
|
||||
|
||||
|
||||
SUBROUTINE par_ptdslv(r,la,wa,nd,n,nall)
|
||||
INTEGER,PARAMETER :: nb = 64
|
||||
INTEGER,INTENT(IN) :: n,nd,nall
|
||||
|
|
|
|||
|
|
@ -546,8 +546,8 @@ module m_fdm_calc
|
|||
ENDDO
|
||||
ENDDO
|
||||
|
||||
CALL par_dfp (fz, fzl, fzu, dfz, hy, 4*xx*yy, zz, yy, 3)
|
||||
CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, xx*yy, zz, yy, 3)
|
||||
CALL pdfp (fz, fzl, fzu, dfz, hy, 4*xx*yy, zz, yy, 3)
|
||||
CALL pd2fp(fzz, fzzl, fzzu, dfzz, hy, xx*yy, zz, yy, 3)
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue