From e659052e26ff8885e5b844405007972f7e615b8e Mon Sep 17 00:00:00 2001 From: ignis Date: Fri, 10 Nov 2017 23:39:52 +0900 Subject: [PATCH] z-dir rhs exchange mpi subroutine change, rk substep advance opt --- m_compact.f90 | 141 +++++++++++++++++++++++------------------------- m_fdm_calc.f90 | 142 ++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 184 insertions(+), 99 deletions(-) diff --git a/m_compact.f90 b/m_compact.f90 index 4e83868..3b789bf 100644 --- a/m_compact.f90 +++ b/m_compact.f90 @@ -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 diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 24b236d..90862e8 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -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