From 3225cd757bad102bd4ac1f94a63b892f39fa1f84 Mon Sep 17 00:00:00 2001 From: ignis Date: Mon, 6 Feb 2017 04:02:33 +0900 Subject: [PATCH] fully mpi working --- m_compact.f90 | 156 ++++++++++++++++++++++++++++++++----------------- m_fdm_calc.f90 | 110 +++++++++++++++------------------- m_mpi.f90 | 3 +- m_openmpi.f90 | 3 +- 4 files changed, 153 insertions(+), 119 deletions(-) diff --git a/m_compact.f90 b/m_compact.f90 index 08fd1fa..eb98435 100644 --- a/m_compact.f90 +++ b/m_compact.f90 @@ -222,7 +222,7 @@ END SUBROUTINE dfp - SUBROUTINE par_dfp(x,dx,h,nd,n,nall,dir) + SUBROUTINE par_dfp(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 @@ -295,7 +295,6 @@ ENDDO - call MPI_BARRIER(MPI_COMM_TASK, mpi_err) IF (dir.eq.1) CALL par_ptdslv(dx,lxf,wxf,nd,n,nall) ! x-direction IF (dir.eq.2) CALL par_ptdslv(dx,lyf,wyf,nd,n,nall) ! x-direction IF (dir.eq.3) CALL par_ptdslv(dx,lzf,wzf,nd,n,nall) ! x-direction @@ -330,17 +329,20 @@ ENDDO END SUBROUTINE ptdslv - SUBROUTINE par_ptdslv(r,l,w,nd,n,nall) + SUBROUTINE par_ptdslv(r,la,wa,nd,n,nall) + INTEGER,PARAMETER :: nb = 128 INTEGER,INTENT(IN) :: n,nd,nall REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r - REAL*8,INTENT(IN),DIMENSION(:) :: l,w + REAL*8,INTENT(IN),DIMENSION(nall) :: la,wa + REAL*8,DIMENSION(n) :: l,w INTEGER i,j INTEGER ii,jj REAL*8, DIMENSION(nd) :: sum INTEGER npart, nbase, nlow, nupp - INTEGER :: pid,np - REAL*8 :: r0, r1, sum0 + INTEGER :: pid, np + REAL*8, DIMENSION(nb) :: r0, r1, sum0 + REAL*8, DIMENSION(nb,2) :: buf pid = myid np = numprocs @@ -350,6 +352,9 @@ nlow = pid * npart + 1 nupp = (pid + 1) * npart + l = la(nlow:nupp) + w = wa(nlow:nupp) + if (npart.lt.4) then ! assertion fail endif @@ -357,28 +362,42 @@ call MPI_BARRIER(MPI_COMM_WORLD, mpi_err) ! first process - if (myid.eq.master) then + if (myid.eq.0) then - DO j=1,nd + DO jj=1,nd,nb + DO j=jj,jj+nb-1 sum(j)=w(1)*r(j,1) r(j,1)=r(j,1)*l(1) + ENDDO - DO i=2,npart + DO i=2,n + DO j=jj,jj+nb-1 r(j,i)=r(j,i)-r(j,i-1) sum(j)=sum(j)+w(i)*r(j,i) r(j,i)=r(j,i)*l(i) + ENDDO ENDDO - CALL MPI_ISEND(r(j,npart), 1, MPI_REAL8, pid+1, 0, MPI_COMM_TASK, mpi_request, mpi_err) - CALL MPI_ISEND(sum(j), 1, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_request, mpi_err) + buf(:,1) = r(jj:jj+nb-1,n) + buf(:,2) = sum(jj:jj+nb-1) + CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err) - CALL MPI_RECV(r0, 1, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err) - CALL MPI_RECV(r1, 1, MPI_REAL8, pid+1, 3, MPI_COMM_TASK, mpi_status, mpi_err) + ENDDO - i=npart - r(j,i)=r(j,i)-l(i)*r0-w(i)*r1 - DO i=npart-1,1,-1 - r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1 + DO jj=1,nd,nb + + CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err) + r0 = buf(:,1) + r1 = buf(:,2) + + i=n + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1) + ENDDO + DO i=n-1,1,-1 + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1) + ENDDO ENDDO ENDDO @@ -386,68 +405,100 @@ ! last process elseif (pid.eq.(np-1)) then - DO j=1,nd + DO jj=1,nd,nb - CALL MPI_RECV(r0, 1, MPI_REAL8, pid-1, 0, MPI_COMM_TASK, mpi_status, mpi_err) - CALL MPI_RECV(sum0, 1, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err) + CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err) + r0 = buf(:,1) + sum0 = buf(:,2) i=1 - r(j,i)=r(j,i)-r0 - sum(j)=sum0+w(i+nbase)*r(j,i) - r(j,i)=r(j,i)*l(i+nbase) + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-r0(j-jj+1) + sum(j)=sum0(j-jj+1)+w(i)*r(j,i) + r(j,i)=r(j,i)*l(i) + ENDDO - DO i=2,npart-1 + DO i=2,n-1 + DO j=jj,jj+nb-1 r(j,i)=r(j,i)-r(j,i-1) - sum(j)=sum(j)+w(i+nbase)*r(j,i) - r(j,i)=r(j,i)*l(i+nbase) + sum(j)=sum(j)+w(i)*r(j,i) + r(j,i)=r(j,i)*l(i) + ENDDO ENDDO - r(j,npart)=l(nupp)*(r(j,npart)-sum(j)) + DO j=jj,jj+nb-1 + r(j,n)=l(n)*(r(j,n)-sum(j)) + ENDDO - r(j,npart-1)=r(j,npart-1)-w(nupp-1)*r(j,npart) + ENDDO - DO i=npart-2,1,-1 - r(j,i)=r(j,i)-l(i+nbase)*r(j,i+1)-w(i+nbase)*r(j,npart) + DO jj=1,nd,nb + + DO j=jj,jj+nb-1 + r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n) + ENDDO + + DO i=n-2,1,-1 + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n) + ENDDO ENDDO - CALL MPI_ISEND(r(j,1), 1, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_request, mpi_err) - CALL MPI_ISEND(r(j,npart), 1, MPI_REAL8, pid-1, 3, MPI_COMM_TASK, mpi_request, mpi_err) + buf(:,1) = r(jj:jj+nb-1,1) + buf(:,2) = r(jj:jj+nb-1,n) + CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err) ENDDO ! intermediate process else - DO j=1,nd + DO jj=1,nd,nb - CALL MPI_RECV(r0, 1, MPI_REAL8, pid-1, 0, MPI_COMM_TASK, mpi_status, mpi_err) - CALL MPI_RECV(sum0, 1, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err) + CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err) + r0 = buf(:,1) + sum0 = buf(:,2) i=1 - r(j,i)=r(j,i)-r0 - sum(j)=sum0+w(i+nbase)*r(j,i) - r(j,i)=r(j,i)*l(i+nbase) + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-r0(j-jj+1) + sum(j)=sum0(j-jj+1)+w(i)*r(j,i) + r(j,i)=r(j,i)*l(i) + ENDDO - DO i=2,npart + DO i=2,n + DO j=jj,jj+nb-1 r(j,i)=r(j,i)-r(j,i-1) - sum(j)=sum(j)+w(i+nbase)*r(j,i) - r(j,i)=r(j,i)*l(i+nbase) + sum(j)=sum(j)+w(i)*r(j,i) + r(j,i)=r(j,i)*l(i) + ENDDO ENDDO - CALL MPI_ISEND(r(j,npart), 1, MPI_REAL8, pid+1, 0, MPI_COMM_TASK, mpi_request, mpi_err) - CALL MPI_ISEND(sum(j), 1, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_request, mpi_err) + buf(:,1) = r(jj:jj+nb-1,n) + buf(:,2) = sum(jj:jj+nb-1) + CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err) - CALL MPI_RECV(r0, 1, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err) - CALL MPI_RECV(r1, 1, MPI_REAL8, pid+1, 3, MPI_COMM_TASK, mpi_status, mpi_err) + ENDDO - i=npart - r(j,i)=r(j,i)-l(i+nbase)*r0-w(i+nbase)*r1 - DO i=npart-1,1,-1 - r(j,i)=r(j,i)-l(i+nbase)*r(j,i+1)-w(i+nbase)*r1 + DO jj=1,nd,nb + + CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err) + r0 = buf(:,1) + r1 = buf(:,2) + + i=n + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1) + ENDDO + DO i=n-1,1,-1 + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1) + ENDDO ENDDO - CALL MPI_ISEND(r(j,1), 1, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_request, mpi_err) - CALL MPI_ISEND(r1, 1, MPI_REAL8, pid-1, 3, MPI_COMM_TASK, mpi_request, mpi_err) + buf(:,1) = r(jj:jj+nb-1,1) + !buf(:,2) = r1 + CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err) ENDDO @@ -518,7 +569,7 @@ IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction END SUBROUTINE d2fp - SUBROUTINE par_d2fp(x,dx,h,nd,n,nall,dir) + SUBROUTINE par_d2fp(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 @@ -626,7 +677,6 @@ ENDDO ENDDO - call MPI_BARRIER(MPI_COMM_TASK, mpi_err) IF (dir.eq.1) CALL par_ptdslv(dx,lxs,wxs,nd,n,nall) ! x-direction IF (dir.eq.2) CALL par_ptdslv(dx,lys,wys,nd,n,nall) ! y-direction IF (dir.eq.3) CALL par_ptdslv(dx,lzs,wzs,nd,n,nall) ! z-direction diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 75d3682..375ca44 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -10,6 +10,8 @@ module m_fdm_calc real*8, dimension(:,:,:), allocatable :: u_,v_,w_ real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf + real*8, dimension(:,:), allocatable :: fz, dfz, fzz, dfzz + real*8, dimension(:,:), allocatable :: fzl, fzu, fzzl, fzzu real*8 :: in_yr,out_yr,refwr,minf integer :: fullsavenum !,svfx,svfy @@ -23,11 +25,6 @@ module m_fdm_calc real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon real*8, dimension(3) :: velmax, velmin, velmax1, velmin1 - integer :: vtype1, vtype2 - integer, dimension(128) :: scnt1, sdisp1, stype1, rcnt1, rdisp1, rtype1 - integer, dimension(128) :: scnt2, sdisp2, stype2, rcnt2, rdisp2, rtype2 - - !=========================================================================== !=========================================================================== @@ -241,40 +238,23 @@ module m_fdm_calc allocate(y2(2,nx,ny,nz)) allocate(yf(2,nx,ny,nz)) + allocate(fz(4*nx*ny,nz)) + allocate(dfz(4*nx*ny,nz)) + allocate(fzz(nx*ny,nz)) + allocate(dfzz(nx*ny,nz)) + + allocate(fzu(4*nx*ny,2)) + allocate(fzl(4*nx*ny,2)) + + allocate(fzzu(nx*ny,2)) + allocate(fzzl(nx*ny,2)) + y1=0.0 y2=0.0 yf=0.0 CALL ludcmp(nx,ny,nz_all,1,0,0) - CALL MPI_TYPE_VECTOR (nz, nx*nz, nx*ny, MPI_REAL8, vtype1, mpi_err) - CALL MPI_TYPE_COMMIT (vtype1, mpi_err) - - CALL MPI_TYPE_VECTOR (nz, 2*nx*nz, 2*nx*ny, MPI_REAL8, vtype2, mpi_err) - CALL MPI_TYPE_COMMIT (vtype2, mpi_err) - - - do i = 1, (ny/nz) - scnt1(i) = 1 - rcnt1(i) = nx * nz * nz - - sdisp1(i) = (i-1) * nx * nz * 8 - rdisp1(i) = (i-1) * nx * nz * nz * 8 - - stype1(i) = vtype1 - rtype1(i) = MPI_REAL8 - - scnt2(i) = 1 - rcnt2(i) = 2 * nx * nz * nz - - sdisp2(i) = (i-1) * 2 * nx * nz * 8 - rdisp2(i) = (i-1) * 2 * nx * nz * nz * 8 - - stype2(i) = vtype2 - rtype2(i) = MPI_REAL8 - - end do - refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! Kwon @@ -295,8 +275,8 @@ module m_fdm_calc in_yr=yy(1) ! inlet_Yr out_yr=yy(nx) ! outlet_Yr - do i=1,ny - do j=1,nz + do j=1,nz + do i=1,ny do ii=1,nx y1(1,ii,i,j)=1. ! rho initializing y1(2,ii,i,j)=yy(ii) ! Yr initializing @@ -332,8 +312,8 @@ module m_fdm_calc CALL MPI_ALLREDUCE(sumc,sumc1,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) fl_location=(hx*(nx-1.))*(1.-(sumc1/(nx*ny*nz_all))) - sumc1=sumc*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.)) - CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) + sumc=sumc1*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.)) + !CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) oldsumc=sumc if (myid.eq.0) write(*,634) fl_location/(REAL(nx-1)*hx)*100. @@ -487,16 +467,17 @@ module m_fdm_calc integer :: i,j,k,xx,yy,zz,ii integer :: n + integer :: idx1, idx2 real*8 :: wrate,yr,yp real*8 :: r1_(2,xx,yy,zz),f_(2,xx,yy,zz) real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: ux(4,xx),dux(4,xx),d2ux(xx) real*8 :: uy(4,yy),duy(4,yy),d2uy(yy) real*8 :: uz(4,zz),duz(4,zz),d2uz(zz) - !real*8 :: uzz(4*xx,yy),duzz(4*xx,yy),d2uzz(xx,yy) real*8 :: uux(xx) real*8 :: uuy(yy) real*8 :: uuz(zz) + real*8 :: y ! reaction source term @@ -504,12 +485,12 @@ module m_fdm_calc DO j=1,yy DO i=1,xx - ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y + y=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y - wrate=pre*ux(2,i)*exp(-ac/(1.+bc*(1.-ux(2,i)))) !wrate - IF ((1.-ux(2,i)).le.c_ref) THEN + wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate + IF ((1.-y).le.c_ref) THEN wrate=min_wr - IF ((1.-ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(2,i)-c_ref))+ & + IF ((1.-y).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-y-c_ref))+ & min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) ENDIF @@ -544,42 +525,43 @@ module m_fdm_calc ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz ! = -( d(rho*w*Yr)/dz ) ! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz ) - f_(2,i,j,k)= f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv. + f_(2,i,j,k) = f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv. ENDDO ENDDO ENDDO else - DO j=1,yy - !DO k=1,zz - ! DO i=1,xx - ! uzz(0*xx+i,k)=r1_(1,i,j,k) ! 1:rho - ! uzz(1*xx+i,k)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y - ! uzz(2*xx+i,k)=uzz(0*xx+i,k)*ww_(i,j,k) ! 3:rho*w - ! uzz(3*xx+i,k)=uz(2*xx+i,k)*uz(1*xx+i,k) ! 4:rho*w*Y - ! ENDDO - !ENDDO - DO i=1,xx - DO k=1,zz - uz(1,k)=r1_(1,i,j,k) ! 1:rho - uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y - uz(3,k)=uz(1,k)*ww_(i,j,k) ! 3:rho*w - uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*Y - uuz (k)=uz(2,k) + DO k=1,zz + DO j=1,yy + DO i=1,xx + idx2 = xx*(j-1)+i + idx1 = (idx2-1)*4 + fz(idx1+1,k) = r1_(1,i,j,k) ! 1:rho + fz(idx1+2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y + fz(idx1+3,k) = r1_(1,i,j,k)*ww_(i,j,k) ! 3:rho*w + fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k) ! 4:rho*w*Y + fzz(idx2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ENDDO + ENDDO + ENDDO - CALL par_dfp (uz(1:4,:),duz(1:4,:),hy,4,zz,yy,3) - CALL par_d2fp(uuz(:), d2uz(:), hy,1,zz,yy,3) + 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) - DO k=1,zz + DO k=1,zz + DO j=1,yy + DO i=1,xx + idx2 = xx*(j-1)+i + idx1 = (idx2-1)*4 ! -( d(rho*w)/dz ) - f_(1,i,j,k) = f_(1,i,j,k) - duz(3,k) ! continuity + f_(1,i,j,k) = f_(1,i,j,k) - dfz(idx1+3,k) ! continuity ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz ! = -( d(rho*w*Yr)/dz ) ! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz ) - f_(2,i,j,k)= f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv. + f_(2,i,j,k) = f_(2,i,j,k) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv. + ENDDO ENDDO ENDDO diff --git a/m_mpi.f90 b/m_mpi.f90 index 3cc221a..b3117c6 100644 --- a/m_mpi.f90 +++ b/m_mpi.f90 @@ -64,7 +64,8 @@ contains ! initializing MPI environment - call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err) + !call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err) + call MPI_INIT(mpi_err) call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err) call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err) diff --git a/m_openmpi.f90 b/m_openmpi.f90 index 73d297a..97ed064 100644 --- a/m_openmpi.f90 +++ b/m_openmpi.f90 @@ -64,7 +64,8 @@ contains ! initializing MPI environment - call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err) + !call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err) + call MPI_INIT(mpi_err) call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err) call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err)