diff --git a/Makefile b/Makefile index 0572750..4a40277 100644 --- a/Makefile +++ b/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 diff --git a/m_compact.f90 b/m_compact.f90 index fd175d6..08fd1fa 100644 --- a/m_compact.f90 +++ b/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 diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index c3ded7d..75d3682 100644 --- a/m_fdm_calc.f90 +++ b/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 )