z-dir rhs exchange mpi subroutine change, rk substep advance opt

This commit is contained in:
ignis 2017-11-10 23:39:52 +09:00
parent b426be89ad
commit e659052e26
2 changed files with 184 additions and 99 deletions

View file

@ -203,8 +203,8 @@
REAL*8 :: r1,r2,h1
integer (kind=MPI_INTEGER_KIND) :: idx
integer (kind=MPI_INTEGER_KIND),dimension(2) :: requests
logical :: flagLow, flagUpp
integer (kind=MPI_INTEGER_KIND),dimension(4) :: requests
logical :: recvLow, recvUpp
h1=1./h
r1=7./3.
@ -213,62 +213,62 @@
!if (myid.eq.0) write(*,*) "parallel dfp commuication"
if (myid.eq.master) then
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_IRECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_IRECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(4), mpi_err)
elseif (myid.eq.numprocs-1) then
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_IRECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_IRECV (xu, 2*nd, MPI_REAL8, 0, 0, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, 0, 0, &
MPI_COMM_TASK, requests(4), mpi_err)
else
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_IRECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_IRECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(4), mpi_err)
endif
call MPI_STARTALL (4, requests, mpi_err)
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
flagLow = .true.
flagUpp = .true.
recvLow = .true.
recvUpp = .true.
do while (flagLow.or.flagUpp)
do while (recvLow.or.recvUpp)
call MPI_WAITANY (2, requests, idx, mpi_status, mpi_err)
!call MPI_WAIT (request1, mpi_status, mpi_err)
call MPI_WAITANY (2, requests(3:), idx, mpi_status, mpi_err)
if (idx.eq.1) then
DO j=1,nd
@ -279,11 +279,9 @@
dx(j,2) =(r1*( x(j,3)- x(j,1)) +r2*( x(j,4)-xl(j,2))) *h1
ENDDO
flagLow = .false.
recvLow = .false.
endif
!call MPI_WAIT (request2, mpi_status, mpi_err)
if (idx.eq.2) then
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
@ -293,7 +291,7 @@
dx(j,n) =(r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
ENDDO
flagUpp = .false.
recvUpp = .false.
endif
enddo
@ -315,8 +313,8 @@
REAL*8,DIMENSION(nd,2) :: xu, xl
integer (kind=MPI_INTEGER_KIND) :: idx
integer (kind=MPI_INTEGER_KIND),dimension(2) :: requests
logical :: flagLow, flagUpp
integer (kind=MPI_INTEGER_KIND),dimension(4) :: requests
logical :: recvLow, recvUpp
h2=1./(h*h)
r1=6.
@ -325,48 +323,50 @@
!if (myid.eq.0) write(*,*) "parallel d2fp commuication"
if (myid.eq.master) then
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_IRECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_IRECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(4), mpi_err)
elseif (myid.eq.numprocs-1) then
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_IRECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_IRECV (xu, 2*nd, MPI_REAL8, 0, 0, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, 0, 0, &
MPI_COMM_TASK, requests(4), mpi_err)
else
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_IRECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_IRECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, requests(4), mpi_err)
endif
call MPI_STARTALL (4, requests, mpi_err)
DO i=3,n-2
DO j=1,nd
@ -377,15 +377,12 @@
ENDDO
ENDDO
recvLow = .true.
recvUpp = .true.
flagLow = .true.
flagUpp = .true.
do while (recvLow.or.recvUpp)
do while (flagLow.or.flagUpp)
call MPI_WAITANY (2, requests, idx, mpi_status, mpi_err)
!call MPI_WAIT (request1, mpi_status, mpi_err)
call MPI_WAITANY (2, requests(3:), idx, mpi_status, mpi_err)
if (idx.eq.1) then
@ -405,11 +402,9 @@
ENDDO
flagLow = .false.
recvLow = .false.
endif
!call MPI_WAIT (request2, mpi_status, mpi_err)
if (idx.eq.2) then
DO j=1,nd
@ -428,7 +423,7 @@
ENDDO
flagUpp = .false.
recvUpp = .false.
endif
enddo

View file

@ -712,11 +712,11 @@ module m_fdm_calc
real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=2; CALL substep(yy1,yy2,yy1,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=3; CALL substep(yy2,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=4; CALL substep(yy1,yy2,yy1,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=5; CALL substep(yy2,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
CALL rotarysubstep(1, xx, yy, zz, uu_, vv_, ww_, yy1, yy1, yy2)
CALL rotarysubstep(2, xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
CALL rotarysubstep(3, xx, yy, zz, uu_, vv_, ww_, yy2, yyf, yy1)
CALL rotarysubstep(4, xx, yy, zz, uu_, vv_, ww_, yyf, yy1, yy2)
CALL lastsubstep (xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf)
return
END SUBROUTINE RK4
@ -735,6 +735,96 @@ module m_fdm_calc
return
END SUBROUTINE solve
subroutine rotarysubstep(istage,xx,yy,zz,uu_,vv_,ww_,r0,r1,r2)
implicit none
integer :: xx,yy,zz,istage
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
real*8 :: r0(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq)
real*8 :: a(5),b(5)
real*8 :: at, bt
integer :: i,j,k,l
a(1)= 970286171893./4311952581923.
a(2)= 6584761158862./12103376702013.
a(3)= 2251764453980./15575788980749.
a(4)= 26877169314380./34165994151039.
a(5)=0.
b(1)= 1153189308089./22510343858157.
b(2)= 1772645290293./4653164025191.
b(3)= -1672844663538./4480602732383.
b(4)= 2114624349019./3568978502595.
b(5)= 5198255086312./14908931495163.
CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_)
at=a(istage)*fdmdt
bt=(b(istage) - a(istage))*fdmdt
DO l = 1,neq
DO k = 1,zz
DO j = 1,yy
DO i = 1,xx
r1(i,j,k,l) = r1(i,j,k,l) + at*r2(i,j,k,l)
r2(i,j,k,l) = r1(i,j,k,l) + bt*r2(i,j,k,l)
ENDDO
ENDDO
ENDDO
ENDDO
return
END SUBROUTINE rotarysubstep
subroutine lastsubstep(xx,yy,zz,uu_,vv_,ww_,r0,r1,r2)
implicit none
integer :: xx, yy, zz
real*8 :: bt
real*8 :: r0(xx,yy,zz,neq), r1(xx,yy,zz,neq), r2(xx,yy,zz,neq)
real*8 :: uu_(xx,yy,zz), vv_(xx,yy,zz), ww_(xx,yy,zz)
real*8 :: b
integer :: i, j, k
b = 5198255086312./14908931495163.
CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_)
bt = b*fdmdt
r0 = r1 + bt*r2
!==========rho=1 treatment
r0(:,:,:,2) = r0(:,:,:,2)/r0(:,:,:,1)
r0(:,:,:,1) = 1.
DO k = 1,zz
DO j = 1,yy
DO i = 1,xx
!==========Max Yr=1 treatment
r0(i,j,k,2)=MIN(in_yr,r0(i,j,k,2))
!==========Min Yr=0 treatment
! r0(i,j,k,2)=MAX(out_yr,r0(i,j,k,2))
ENDDO
ENDDO
ENDDO
return
END SUBROUTINE lastsubstep
subroutine substep(ri,r1,r2,f,xx,yy,zz,istage,uu_,vv_,ww_)
implicit none
@ -752,6 +842,7 @@ module m_fdm_calc
a(3)= 2251764453980./15575788980749.
a(4)= 26877169314380./34165994151039.
a(5)=0.
b(1)= 1153189308089./22510343858157.
b(2)= 1772645290293./4653164025191.
b(3)= -1672844663538./4480602732383.
@ -761,38 +852,37 @@ module m_fdm_calc
CALL fns(ri,f,xx,yy,zz,uu_,vv_,ww_)
IF(istage<5) THEN
at=a(istage)*fdmdt
bt=(b(istage)-a(istage))*fdmdt
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,neq
r1(i,j,k,nv)=r1(i,j,k,nv)+at*f(i,j,k,nv)
r2(i,j,k,nv)=r1(i,j,k,nv)+bt*f(i,j,k,nv)
ENDDO
ENDDO
ENDDO
ENDDO
ELSE
bt=b(istage)*fdmdt
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,neq
r1(i,j,k,nv)=r1(i,j,k,nv)+bt*f(i,j,k,nv)
ENDDO
r2 = r1 + bt*f
r1 = r1 + at*f
ELSE
bt=b(istage)*fdmdt
r1 = r1 + bt*f
!==========rho=1 treatment
r1(i,j,k,2)=r1(i,j,k,2)/r1(i,j,k,1)
r1(i,j,k,1)=1.
r1(:,:,:,2) = r1(:,:,:,2)/r1(:,:,:,1)
r1(:,:,:,1) = 1.
DO k = 1,zz
DO j = 1,yy
DO i = 1,xx
!==========Max Yr=1 treatment
r1(i,j,k,2)=MIN(in_yr,r1(i,j,k,2))
!==========Min Yr=0 treatment
! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))
! r1(i,j,k,2)=MAX(out_yr,r1(i,j,k,2))
ENDDO
ENDDO
ENDDO
ENDIF
return