research/dns-hit3d-fdm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

fully mpi working

Browse source
This commit is contained in:
ignis 2017-02-06 04:02:33 +09:00
parent 6c75483f1e
commit 3225cd757b
4 changed files with 153 additions and 119 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

156
m_compact.f90
View file

@ -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

110
m_fdm_calc.f90
View file

@ -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

3
m_mpi.f90
View file

@ -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)

3
m_openmpi.f90
View file

@ -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)