pipeline parallel compact scheme
This commit is contained in:
ignis
parent
383b2c5aa7
commit
6c75483f1e
3 changed files with 417 additions and 122 deletions
62
Makefile
62
Makefile
|
|
@ -4,7 +4,6 @@
|
|||
|
||||
MPIF90 = blah
|
||||
|
||||
ifeq ($(OMPI), 1)
|
||||
FFTW_HOME = ../libs/fftw-3.2.2/lib
|
||||
OPENMPI_HOME = /home/comb/ignis/openmpi/
|
||||
MPIF90 = $(OPENMPI_HOME)/bin/mpif90
|
||||
|
|
@ -13,67 +12,6 @@ ifeq ($(OMPI), 1)
|
|||
FCFLAGS_F77 = -Mextend
|
||||
|
||||
MODULES += m_openmpi.o
|
||||
else
|
||||
FFTW_HOME = ../libs/fftw-3.2.2/lib
|
||||
MPIF90 = mpif90
|
||||
FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
|
||||
LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
|
||||
FCFLAGS_F77 = -Mextend
|
||||
|
||||
MODULES += m_mpi.o
|
||||
endif
|
||||
|
||||
# Franklin cluster, NERSC
|
||||
ifeq ($(NERSC_HOST), franklin)
|
||||
MPIF90 = ftn
|
||||
|
||||
FCFLAGS = -target=linux -i8 -r8 -O4 -c
|
||||
LDFLAGS = -target=linux -i8 -r8 -O4 -lfftw3
|
||||
FCFLAGS_F77 = -Mextend
|
||||
endif
|
||||
|
||||
# Hopper cluster, NERSC
|
||||
ifeq ($(NERSC_HOST), hopper)
|
||||
MPIF90 = ftn
|
||||
|
||||
FCFLAGS = -target=linux -i8 -r8 -O4 -c
|
||||
LDFLAGS = -target=linux -i8 -r8 -O4 -lfftw3
|
||||
FCFLAGS_F77 = -Mextend
|
||||
endif
|
||||
|
||||
# Yellowrail cluster, LANL
|
||||
ifeq ($(HOSTNAME), yr-fe1.lanl.gov)
|
||||
MPIF90 = mpif90
|
||||
|
||||
FCFLAGS = -i8 -r8 -O4 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_INCLUDE)
|
||||
LDFLAGS = -i8 -r8 -O4 -lfftw3 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib
|
||||
endif
|
||||
|
||||
# Coyote cluster, LANL
|
||||
ifeq ($(HOSTNAME), cy-c2.lanl.gov)
|
||||
MPIF90 = mpif90
|
||||
FCFLAGS = -i8 -r8 -O4 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_INCLUDE)
|
||||
LDFLAGS = -i8 -r8 -O4 -lfftw3 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib
|
||||
endif
|
||||
|
||||
# Linux Cluster, CMU
|
||||
ifeq ($(HOSTNAME), karman.me.cmu.edu)
|
||||
FFTW_HOME = ../libs/fftw-3.2.2/lib
|
||||
MPIF90 = mpif90
|
||||
FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
|
||||
LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
|
||||
FCFLAGS_F77 = -Mextend
|
||||
|
||||
endif
|
||||
|
||||
# sparrow.stanford.edu, Mac OS X box
|
||||
ifeq ($(HOSTNAME), sparrow.stanford.edu)
|
||||
MPIF90 = mpif90
|
||||
FCFLAGS = -fdefault-real-8 -fdefault-integer-8 -finit-integer=0 -finit-real=zero -c
|
||||
FCFLAGS_F77 = -ffixed-form -ffixed-line-length-none
|
||||
FCFLAGS_F90 = -ffree-form -ffree-line-length-none
|
||||
LDFLAGS = -fdefault-real-8 -fdefault-integer-8 -finit-integer=0 -finit-real=zero -lmpi -lmpi -lfftw3 -lm
|
||||
endif
|
||||
|
||||
# Program name
|
||||
PROG = hit3d.x
|
||||
|
|
|
|||
324
m_compact.f90
324
m_compact.f90
|
|
@ -1,4 +1,7 @@
|
|||
MODULE m_compact
|
||||
|
||||
use m_openmpi
|
||||
|
||||
IMPLICIT NONE
|
||||
PRIVATE
|
||||
REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
|
||||
|
|
@ -8,7 +11,7 @@
|
|||
INTEGER :: nxc,nyc,nzc
|
||||
REAL*8, PARAMETER :: ezero = 1.0e-14
|
||||
|
||||
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
|
||||
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
@ -219,6 +222,86 @@
|
|||
|
||||
END SUBROUTINE dfp
|
||||
|
||||
SUBROUTINE par_dfp(x,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 :: r1,r2,h1
|
||||
|
||||
REAL*8,DIMENSION(nd,2) :: xu, xl
|
||||
|
||||
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,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
|
||||
dx(j,n)= (r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
|
||||
dx(j,1)= (r1*( x(j,2)-xl(j,2)) +r2*( x(j,3)-xl(j,1)))*h1
|
||||
dx(j,2)= (r1*( x(j,3)- x(j,1)) +r2*( x(j,4)-xl(j,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
|
||||
|
||||
|
||||
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
|
||||
|
||||
END SUBROUTINE par_dfp
|
||||
|
||||
SUBROUTINE ptdslv(r,n,l,w,nd)
|
||||
INTEGER,INTENT(IN) :: n,nd
|
||||
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
|
|
@ -247,6 +330,131 @@
|
|||
ENDDO
|
||||
END SUBROUTINE ptdslv
|
||||
|
||||
SUBROUTINE par_ptdslv(r,l,w,nd,n,nall)
|
||||
INTEGER,INTENT(IN) :: n,nd,nall
|
||||
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
REAL*8,INTENT(IN),DIMENSION(:) :: 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
|
||||
|
||||
pid = myid
|
||||
np = numprocs
|
||||
|
||||
npart = nall / np
|
||||
nbase = pid * npart
|
||||
nlow = pid * npart + 1
|
||||
nupp = (pid + 1) * npart
|
||||
|
||||
if (npart.lt.4) then
|
||||
! assertion fail
|
||||
endif
|
||||
|
||||
call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)
|
||||
|
||||
! first process
|
||||
if (myid.eq.master) then
|
||||
|
||||
DO j=1,nd
|
||||
sum(j)=w(1)*r(j,1)
|
||||
r(j,1)=r(j,1)*l(1)
|
||||
|
||||
DO i=2,npart
|
||||
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
|
||||
|
||||
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)
|
||||
|
||||
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)
|
||||
|
||||
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
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
|
||||
! last process
|
||||
elseif (pid.eq.(np-1)) then
|
||||
|
||||
DO j=1,nd
|
||||
|
||||
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)
|
||||
|
||||
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 i=2,npart-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)
|
||||
ENDDO
|
||||
|
||||
r(j,npart)=l(nupp)*(r(j,npart)-sum(j))
|
||||
|
||||
r(j,npart-1)=r(j,npart-1)-w(nupp-1)*r(j,npart)
|
||||
|
||||
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)
|
||||
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)
|
||||
|
||||
ENDDO
|
||||
|
||||
! intermediate process
|
||||
else
|
||||
|
||||
DO j=1,nd
|
||||
|
||||
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)
|
||||
|
||||
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 i=2,npart
|
||||
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)
|
||||
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)
|
||||
|
||||
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)
|
||||
|
||||
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
|
||||
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)
|
||||
|
||||
ENDDO
|
||||
|
||||
endif
|
||||
|
||||
END SUBROUTINE par_ptdslv
|
||||
|
||||
SUBROUTINE d2fp(n,h,x,dx,nd,dir)
|
||||
INTEGER,INTENT(IN) :: n,nd,dir
|
||||
REAL*8,INTENT(IN) :: h
|
||||
|
|
@ -310,6 +518,120 @@
|
|||
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)
|
||||
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,n)-2.*x(j,n-1)+x(j,n-2))
|
||||
t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
|
||||
IF ( x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0.
|
||||
IF (xu(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0.
|
||||
dx(j,n-1)=(r1*t1+r2*t2)*h2
|
||||
|
||||
t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
|
||||
t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
|
||||
IF (xu(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0.
|
||||
IF (xu(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0.
|
||||
! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
|
||||
! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
|
||||
dx(j,n)=(r1*t1+r2*t2)*h2
|
||||
|
||||
t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
|
||||
t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
|
||||
IF (x(j,2).eq.x(j,1).and.x(j,1).eq.xl(j,2)) t1=0.
|
||||
IF (x(j,3).eq.x(j,1).and.x(j,1).eq.xl(j,1)) t2=0.
|
||||
! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
|
||||
! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
|
||||
dx(j,1)=(r1*t1+r2*t2)*h2
|
||||
|
||||
t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
|
||||
t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
|
||||
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0.
|
||||
IF (x(j,4).eq.x(j,2).and.x(j,2).eq.xl(j,2)) t2=0.
|
||||
! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
|
||||
! +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
|
||||
dx(j,2)=(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))
|
||||
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
|
||||
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
|
||||
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
|
||||
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
|
||||
dx(j,i)=(r1*t1+r2*t2)*h2
|
||||
|
||||
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
|
||||
END SUBROUTINE par_d2fp
|
||||
|
||||
SUBROUTINE tdslv(r,n,l,nd)
|
||||
INTEGER,INTENT(IN) :: n,nd
|
||||
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
|
||||
|
|
|
|||
153
m_fdm_calc.f90
153
m_fdm_calc.f90
|
|
@ -330,7 +330,8 @@ module m_fdm_calc
|
|||
CALL MPI_ALLREDUCE(sum_wrate1,sum_wrate,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
|
||||
if (myid.eq.0) write(*,633) sum_wrate
|
||||
|
||||
fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz_all)))
|
||||
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)
|
||||
oldsumc=sumc
|
||||
|
|
@ -444,7 +445,8 @@ module m_fdm_calc
|
|||
if (myid.eq.0) write(*,633) sum_wrate
|
||||
633 format (' ** Consumption Speed, Sc = ',f7.4)
|
||||
|
||||
fl_location=(hx*REAL(nx-1))*(1.-(sumc/(REAL(nx*ny*nz_all))))
|
||||
CALL MPI_ALLREDUCE(sumc,sumc1,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
|
||||
fl_location=(hx*REAL(nx-1))*(1.-(sumc1/(REAL(nx*ny*nz_all))))
|
||||
sumc1=sumc*(hx*hy*hy)/(hy*(REAL(ny)-1.)*hy*(REAL(ny)-1.))
|
||||
CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
|
||||
if (myid.eq.0) write(*,634) fl_location/(REAL(nx-1)*hx)*100.
|
||||
|
|
@ -484,69 +486,105 @@ module m_fdm_calc
|
|||
implicit none
|
||||
|
||||
integer :: i,j,k,xx,yy,zz,ii
|
||||
integer :: n
|
||||
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,yy),duz(4,yy),d2uz(yy)
|
||||
real*8 :: Ly(yy,zz),Dy
|
||||
|
||||
real*8 :: r1_t(2,xx,zz,yy),f_t(2,xx,zz,yy)
|
||||
real*8 :: ww_t(xx,zz,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)
|
||||
|
||||
|
||||
!CALL MPI_ALLTOALLW (ww_, scnt1, sdisp1, stype1, ww_t, rcnt1, rdisp1, rtype1, MPI_COMM_TASK, mpi_err)
|
||||
! reaction source term
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
|
||||
!CALL MPI_ALLTOALLW (r1_, scnt2, sdisp2, stype2, r1_t, rcnt2, rdisp2, rtype2, MPI_COMM_TASK, mpi_err)
|
||||
ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
|
||||
|
||||
do i = 1,numprocs
|
||||
CALL MPI_ISEND(ww_(1,zz*(i-1)+1,1), scnt1(i), stype1(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
end do
|
||||
do i = 1,numprocs
|
||||
CALL MPI_RECV(ww_t(1,1,zz*(i-1)+1), rcnt1(i), rtype1(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
end do
|
||||
wrate=pre*ux(2,i)*exp(-ac/(1.+bc*(1.-ux(2,i)))) !wrate
|
||||
IF ((1.-ux(2,i)).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))+ &
|
||||
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
||||
ENDIF
|
||||
|
||||
do i = 1,numprocs
|
||||
CALL MPI_ISEND(r1_(1,1,zz*(i-1)+1,1), scnt2(i), stype2(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
end do
|
||||
do i = 1,numprocs
|
||||
CALL MPI_RECV(r1_t(1,1,1,zz*(i-1)+1), rcnt2(i), rtype2(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
end do
|
||||
f_(1,i,j,k) = 0.0 ! continuity
|
||||
f_(2,i,j,k) = - wrate ! species conservation
|
||||
|
||||
!! z-direction
|
||||
DO j=1,zz
|
||||
DO i=1,xx
|
||||
DO k=1,yy
|
||||
uz(1,k)=r1_t(1,i,j,k) ! 1:rho
|
||||
uz(2,k)=r1_t(2,i,j,k)/r1_t(1,i,j,k) ! 2:Y
|
||||
uz(3,k)=uz(1,k)*ww_t(i,j,k) ! 3:rho*w
|
||||
uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*Y
|
||||
ENDDO
|
||||
|
||||
CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3)
|
||||
CALL d2fp(yy,hy,uz(2,:),d2uz(:),1,3)
|
||||
|
||||
DO k=1,yy
|
||||
! -( d(rho*w)/dz )
|
||||
f_t(1,i,j,k)=-duz(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_t(2,i,j,k)=-duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv.
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
!CALL MPI_ALLTOALLW (f_t, rcnt2, rdisp2, rtype2, f_, scnt2, sdisp2, stype2, MPI_COMM_TASK, mpi_err)
|
||||
|
||||
do i = 1,numprocs
|
||||
CALL MPI_ISEND(f_t(1,1,1,zz*(i-1)+1), rcnt2(i), rtype2(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
|
||||
end do
|
||||
do i = 1,numprocs
|
||||
CALL MPI_RECV(f_(1,1,zz*(i-1)+1,1), scnt2(i), stype2(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
|
||||
end do
|
||||
!! z-direction
|
||||
if (numprocs.eq.1) then
|
||||
|
||||
DO j=1,yy
|
||||
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)
|
||||
ENDDO
|
||||
|
||||
CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3)
|
||||
CALL d2fp(yy,hy,uuz(:),d2uz(:),1,3)
|
||||
|
||||
DO k=1,zz
|
||||
! -( d(rho*w)/dz )
|
||||
f_(1,i,j,k) = f_(1,i,j,k) - duz(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.
|
||||
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)
|
||||
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)
|
||||
|
||||
DO k=1,zz
|
||||
! -( d(rho*w)/dz )
|
||||
f_(1,i,j,k) = f_(1,i,j,k) - duz(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.
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
endif
|
||||
|
||||
! x-direction
|
||||
DO k=1,zz
|
||||
|
|
@ -556,19 +594,13 @@ module m_fdm_calc
|
|||
ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
|
||||
ux(3,i)=ux(1,i)*uu_(i,j,k) ! 3:rho*u
|
||||
ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*Y
|
||||
uux (i)=ux(2,i)
|
||||
ENDDO
|
||||
|
||||
CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
|
||||
CALL d2fnonp(xx,hx,ux(2,:),d2ux(:),1,1)
|
||||
CALL d2fnonp(xx,hx,uux(:),d2ux(:),1,1)
|
||||
|
||||
DO i=1,xx
|
||||
wrate=pre*ux(2,i)*exp(-ac/(1.+bc*(1.-ux(2,i)))) !wrate
|
||||
IF ((1.-ux(2,i)).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))+ &
|
||||
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
||||
ENDIF
|
||||
|
||||
|
||||
! -( d(rho*u)/dx )
|
||||
f_(1,i,j,k) = f_(1,i,j,k) - dux(3,i) ! continuity
|
||||
|
|
@ -576,22 +608,25 @@ module m_fdm_calc
|
|||
! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
|
||||
! = -( d(rho*u*Yr)/dx )
|
||||
! + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
|
||||
f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) - wrate ! species conservation
|
||||
f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) ! species conservation
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
!! y-direction
|
||||
DO k=1,zz
|
||||
DO i=1,xx
|
||||
DO j=1,yy
|
||||
uy(1,j)=r1_(1,i,j,k) ! 1:rho
|
||||
uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
|
||||
uy(3,j)=uy(1,j)*vv_(i,j,k) ! 3:rho*v
|
||||
uy(4,j)=uy(3,j)*uy(2,j) ! 4:rho*v*Y
|
||||
uuy (j)=uy(2,j)
|
||||
ENDDO
|
||||
|
||||
CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
|
||||
CALL d2fp(yy,hy,uy(2,:),d2uy(:),1,2)
|
||||
CALL d2fp(yy,hy,uuy(:),d2uy(:),1,2)
|
||||
|
||||
DO j=1,yy
|
||||
! -( d(rho*v)/dy )
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue