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

Compare commits

merge into: research:mpi_fdm
Branches Tags
research:mpi_fdm
research:two-step
research:more-iterations
research:diffusivity
research:test_branch
research:pipeline-compact
research:ppc
...
pull from: research:ppc
Branches Tags
research:two-step
research:more-iterations
research:diffusivity
research:test_branch
research:pipeline-compact
research:ppc
research:mpi_fdm
Sign in to create a new pull request.

16 commits
mpi_fdm ... ppc

Author SHA1 Message Date
park0d
16ec6d36e9 make parameter of number of scalar equations 2017-09-27 04:38:07 +09:00
park0d
e14bb9dd14 disable gather_4 before restart_write -> u' + U written 2017-09-27 04:25:44 +09:00
park0d
fe6444366c add scalar equation, disable reaction 2017-09-27 04:25:44 +09:00
ignis
dac299c060 Makefile update to use gfortran compiler 2017-08-02 19:39:19 +09:00
ignis
3225cd757b fully mpi working 2017-02-06 12:49:15 +09:00
ignis
6c75483f1e pipeline parallel compact scheme 2017-01-24 15:50:20 +09:00
ignis
383b2c5aa7 init_velocity seperate n_shell range test and e_spec value test 2017-01-24 15:50:20 +09:00
ignis
6f0c1f44fa fdm restart for mpi 2017-01-24 15:50:20 +09:00
ignis
4f60757a9a mpi communication established 2017-01-24 15:50:20 +09:00
ignis
016bd9fb80 gather4 only performs inverse transform 2017-01-24 15:50:20 +09:00
ignis
804698f717 every process call prepare_fdm and fdm_exe, do not work correctly 2017-01-24 15:50:20 +09:00
ignis
4f7bb89461 remove omp 2017-01-24 15:50:20 +09:00
ignis
50bc5073c1 remove old HIT3D save after writing fdm save 2017-01-24 15:50:20 +09:00
ignis
9c48452b0f disable io 2017-01-24 15:50:20 +09:00
ignis
6d7b31c555 write tp_field flush after each x loop 2017-01-24 15:03:34 +09:00
CombustionLab.POSTECH
cfe19ffb17 job_runlimit default to one month 2017-01-24 14:09:18 +09:00
11 changed files with 872 additions and 318 deletions
Show stats Expand all files Collapse all files

73
Makefile
View file

@ -4,76 +4,13 @@
MPIF90 = blah
ifeq ($(OMPI), 1)
FFTW_HOME = ../libs/fftw-3.2.2/lib
OPENMPI_HOME = /home/comb/ignis/openmpi/
MPIF90 = $(OPENMPI_HOME)/bin/mpif90
FCFLAGS = -i8 -r8 -c -openmp $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include -I$(OPENMPI_HOME)/include
LDFLAGS = -i8 -r8 -openmp $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -L$(OPENMPI_HOME)/lib -lfftw3 -lm -shared-intel
FCFLAGS_F77 = -Mextend
FFTW_HOME = ../libs/fftw
MPIF90 = mpif90
FCFLAGS = -fdefault-double-8 -fdefault-integer-8 -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
LDFLAGS = -fdefault-double-8 -fdefault-integer-8 -fdefault-real-8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
FCFLAGS_F77 =
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

9
gather_4.f90
View file

@ -45,18 +45,15 @@ subroutine gather_4
! velocities
call xFFT3d(-1,1)
fname = 'u.'//file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
call gather_tmp4(u_)
u_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
call xFFT3d(-1,2)
fname = 'v.'//file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
call gather_tmp4(v_)
v_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
call xFFT3d(-1,3)
fname = 'w.'//file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3)
call gather_tmp4(w_)
w_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3)
! scalars
if (int_scalars) then

4
init_velocity.f90
View file

@ -266,8 +266,10 @@ subroutine init_velocity
do i = 1,nx+2
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .gt. 0 .and. n_shell .le. kmax .and. e_spec(n_shell) .gt. zip) then
if (n_shell .gt. 0 .and. n_shell .le. kmax ) then
if (e_spec(n_shell) .gt. zip) then
fields(i,j,k,1:3) = fields(i,j,k,1:3) * sqrt(e_spec1(n_shell)/e_spec(n_shell))
endif
else
fields(i,j,k,1:3) = zip
end if

1
io_write_4.f90
View file

@ -65,7 +65,6 @@ subroutine io_write_4
wmax = maxval(w_(1:nx,1:ny,1:nz_all))
wmin = minval(w_(1:nx,1:ny,1:nz_all))
call write_tp_field
call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
endif

374
m_compact.f90
View file

@ -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,85 @@
END SUBROUTINE dfp
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
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
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 +329,183 @@
ENDDO
END SUBROUTINE ptdslv
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(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, DIMENSION(nb) :: r0, r1, sum0
REAL*8, DIMENSION(nb,2) :: buf
pid = myid
np = numprocs
npart = nall / np
nbase = pid * npart
nlow = pid * npart + 1
nupp = (pid + 1) * npart
l = la(nlow:nupp)
w = wa(nlow:nupp)
if (npart.lt.4) then
! assertion fail
endif
call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)
! first process
if (myid.eq.0) then
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,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
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)
ENDDO
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
! last process
elseif (pid.eq.(np-1)) then
DO jj=1,nd,nb
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
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,n-1
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
DO j=jj,jj+nb-1
r(j,n)=l(n)*(r(j,n)-sum(j))
ENDDO
ENDDO
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
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 jj=1,nd,nb
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
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,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
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)
ENDDO
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
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
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 +569,119 @@
IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction
END SUBROUTINE d2fp
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
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
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

677
m_fdm_calc.f90
View file

@ -1,24 +1,34 @@
module m_fdm_calc
use m_parameters
use m_io
use m_compact
implicit none
!variables
integer, parameter :: nscalar = 2
integer, parameter :: neq = 1 + nscalar
integer, parameter :: nd1 = 2 + 2 * nscalar
integer, parameter :: nd2 = nscalar
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
integer :: fdmcyc,fdmsavecount
real*8 :: fdmtime,fdmdt
real*8 :: sumc,oldsumc,time_int,sum_wrate
real*8 :: sumc1,sum_wrate1
real*8 :: visdt, convdt
integer :: fdmstep
real*8 :: coe,tt1,tt2,tt3,tt4
real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon
real*8, dimension(3) :: velmax, velmin, velmax1, velmin1
!===========================================================================
!===========================================================================
@ -26,134 +36,28 @@ module m_fdm_calc
contains
subroutine fdm_restart_write
write(*,*) '======================================================='
write(*,*) 'Full results are being written',fullsavenum
OPEN (fullsavenum,form='unformatted',status='unknown')
write (fullsavenum) fdmtime,nx,ny,nz_all,oldsumc,time_int
write (fullsavenum) fdmcyc,DT,dummyu_
write (fullsavenum) dt_fdmsave,dt_fullsave
write (fullsavenum) t_fdmsave,t_fullsave
write (fullsavenum) in_yr,out_yr
write (fullsavenum) u_+dummyu_,v_,w_,y1
CLOSE (fullsavenum)
write(*,*) '======================================================='
fullsavenum=fullsavenum+1
end subroutine fdm_restart_write
subroutine fdm_restart_read
integer(kind=MPI_INTEGER_KIND) :: fh
integer(kind=MPI_OFFSET_KIND) :: offset
subroutine prepare_fdm
integer :: header_length
implicit none
write(fname,"('fort.',i4)") restartnum
integer :: i,j,ii,k
real*8 :: fl_location,wrate,yr,c
real*8 :: tar_sumc
allocate(u_(nx,ny,nz_all))
allocate(v_(nx,ny,nz_all))
allocate(w_(nx,ny,nz_all))
u_=0.0
v_=0.0
w_=0.0
if (fdm_sw .eq. 0) then
return
endif
! DQ initializing
fdmsavecount=1 !FDM save count
sum_wrate=0.
sumc=0. !for adjusting mean velocity
oldsumc=0.
time_int=0.
visdt=9999.
convdt=9999.
fdmcyc=0
fdmtime=0.
fullsavenum=1000 !full save file
allocate(y1(2,nx,ny,nz_all))
allocate(y2(2,nx,ny,nz_all))
allocate(yf(2,nx,ny,nz_all))
CALL ludcmp(nx,ny,nz_all,1,0,0)
refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! Kwon
!FDM normal start==================================
if (restartnum==0) then
OPEN(305,FILE='sfield.bin',form='unformatted',status='unknown')
DO i=1,nx
READ (305) y1(2,i,1,1) ! Yr
ENDDO
CLOSE (305)
in_yr=y1(2,1,1,1) ! inlet_Yr
out_yr=y1(2,nx,1,1) ! outlet_Yr
do i=1,ny
do j=1,nz_all
do ii=1,nx
y1(1,ii,i,j)=1. ! rho initializing
y1(2,ii,i,j)=y1(2,ii,1,1) ! Yr initializing
enddo
enddo
enddo
! flame location setup
sum_wrate=0.; sumc=0.
DO ii=1,nz_all
DO j=1,ny
DO i=1,nx
yr=y1(2,i,j,ii)/y1(1,i,j,ii)
c=1.-yr
IF (c.lt.0.) c=0.
wrate=pre*yr*exp(-ac/(1.+bc*c))
! cold boundary difficulty treatment
IF(c.le.c_ref) THEN
wrate=min_wr
IF(c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
sum_wrate=sum_wrate+wrate*(hx*hy*hy)
sumc=sumc+(1.-yr)
ENDDO
ENDDO
ENDDO
sum_wrate=sum_wrate/(hy*hy*ny*nz_all)
write(*,633) sum_wrate
fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz_all)))
sumc=sumc*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.))
oldsumc=sumc
write(*,634) fl_location/(REAL(nx-1)*hx)*100.
633 format (' ** Consumption Speed, Sc = ',f7.4)
634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
! Restart setup ================================================================
else
if (myid.eq.master) then
write (6,*) '********************************************************'
write (6,*) ' FDM field is being initialized for restarting.'
write (6,457) restartnum
457 format(' Restart file number : ',i6)
OPEN (restartnum,form='unformatted',status='unknown')
OPEN (restartnum,file=fname,form='unformatted',status='unknown')
read (restartnum) fdmtime,tt1,tt2,tt3,oldsumc,time_int
read (restartnum) fdmcyc,tt4,dummyu_
read (restartnum) fdmcyc,fdmdt,dummyu_
read (restartnum) tt1,tt2 !dt_fdmsave,dt_fullsave
read (restartnum) t_fdmsave,t_fullsave
read (restartnum) in_yr,out_yr
read (restartnum) u_,v_,w_,y1
CLOSE (restartnum)
if(fdmtime.ge.t_fdmsave) t_fdmsave=t_fdmsave+dt_fdmsave
if(fdmtime.ge.t_fullsave) t_fullsave=t_fullsave+dt_fullsave
@ -176,6 +80,264 @@ module m_fdm_calc
455 format(' !! Spectral Time : ',f10.5,' / FDM Time : ',f10.5)
endif
end if
count = 1
call MPI_BCAST(fdmtime, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(tt1, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(tt2, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(tt3, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(oldsumc, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(time_int, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(fdmcyc, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(fdmdt, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(dummyu_, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
!call MPI_BCAST(dt_fdmsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
!call MPI_BCAST(dt_fullsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
!call MPI_BCAST(t_fdmsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
!call MPI_BCAST(t_fullsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(in_yr, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
call MPI_BCAST(out_yr, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
! opening the file
call MPI_INFO_CREATE(mpi_info, mpi_err)
call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_RDONLY,mpi_info,fh,mpi_err)
header_length = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8)
count = nx * ny * nz
offset = header_length + 4 + 0*nx*ny*nz_all*8 + myid*nx*ny*nz*8
call MPI_FILE_READ_AT_ALL(fh, offset, u_, count, MPI_REAL8, mpi_status, mpi_err)
offset = header_length + 4 + 1*nx*ny*nz_all*8 + myid*nx*ny*nz*8
call MPI_FILE_READ_AT_ALL(fh, offset, v_, count, MPI_REAL8, mpi_status, mpi_err)
offset = header_length + 4 + 2*nx*ny*nz_all*8 + myid*nx*ny*nz*8
call MPI_FILE_READ_AT_ALL(fh, offset, w_, count, MPI_REAL8, mpi_status, mpi_err)
count = neq * nx * ny * nz
offset = header_length + 4 + 3*nx*ny*nz_all*8 + myid*neq*nx*ny*nz*8
call MPI_FILE_READ_AT_ALL(fh, offset, y1, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_CLOSE(fh, mpi_err)
call MPI_INFO_FREE(mpi_info, mpi_err)
end subroutine fdm_restart_read
subroutine fdm_restart_write
integer(kind=MPI_INTEGER_KIND) :: fh
integer(kind=MPI_OFFSET_KIND) :: offset
integer :: header_length
! --------------- writing process ------------------
if (myid.eq.0) write(*,*) '======================================================='
if (myid.eq.0) write(*,*) 'Full results are being written',fullsavenum
!fname = 'fort.'//savenum
write(fname,"('fort.',i4)") fullsavenum
! opening the file
call MPI_INFO_CREATE(mpi_info, mpi_err)
call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_WRONLY+MPI_MODE_CREATE,mpi_info,fh,mpi_err)
! the master node writes the header with parameters
if (myid.eq.0) then
count = 1
call MPI_FILE_WRITE(fh, 8*6, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, fdmtime, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, nx, count, MPI_INTEGER8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, ny, count, MPI_INTEGER8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, nz_all, count, MPI_INTEGER8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, oldsumc, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, time_int, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*6, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*3, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, fdmcyc, count, MPI_INTEGER8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, DT, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, dummyu_, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*3, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, dt_fdmsave, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, dt_fullsave, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, t_fdmsave, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, t_fullsave, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, in_yr, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, out_yr, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, 8*nx*ny*nz_all*5, count, MPI_INTEGER4, mpi_status, mpi_err)
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 5*nx*ny*nz_all*8
call MPI_FILE_WRITE_AT(fh, offset, 8*nx*ny*nz_all*5, count, MPI_INTEGER4, mpi_status, mpi_err)
end if
header_length = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8)
count = nx * ny * nz
offset = header_length + 4 + 0*nx*ny*nz_all*8 + myid*nx*ny*nz*8
call MPI_FILE_WRITE_AT(fh, offset, u_, count, MPI_REAL8, mpi_status, mpi_err)
offset = header_length + 4 + 1*nx*ny*nz_all*8 + myid*nx*ny*nz*8
call MPI_FILE_WRITE_AT(fh, offset, v_, count, MPI_REAL8, mpi_status, mpi_err)
offset = header_length + 4 + 2*nx*ny*nz_all*8 + myid*nx*ny*nz*8
call MPI_FILE_WRITE_AT(fh, offset, w_, count, MPI_REAL8, mpi_status, mpi_err)
count = neq * nx * ny * nz
offset = header_length + 4 + 3*nx*ny*nz_all*8 + myid*neq*nx*ny*nz*8
call MPI_FILE_WRITE_AT(fh, offset, y1, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_CLOSE(fh, mpi_err)
call MPI_INFO_FREE(mpi_info, mpi_err)
if (myid.eq.0) write(*,*) '======================================================='
fullsavenum=fullsavenum+1
end subroutine fdm_restart_write
subroutine prepare_fdm
implicit none
integer :: i,j,ii,k
real*8 :: fl_location,wrate,yr,c
real*8 :: tar_sumc
real*8 :: yy(nx)
if (fdm_sw .eq. 0) then
return
endif
allocate(u_(nx,ny,nz))
allocate(v_(nx,ny,nz))
allocate(w_(nx,ny,nz))
u_=0.0
v_=0.0
w_=0.0
! DQ initializing
fdmsavecount=1 !FDM save count
sum_wrate=0.
sum_wrate1=0.
sumc=0. !for adjusting mean velocity
sumc1=0. !for adjusting mean velocity
oldsumc=0.
time_int=0.
visdt=9999.
convdt=9999.
fdmcyc=0
fdmtime=0.
fullsavenum=1000 !full save file
allocate(y1(neq,nx,ny,nz))
allocate(y2(neq,nx,ny,nz))
allocate(yf(neq,nx,ny,nz))
allocate( fz (nd1*nx*ny,nz))
allocate(dfz (nd1*nx*ny,nz))
allocate( fzz(nd2*nx*ny,nz))
allocate(dfzz(nd2*nx*ny,nz))
allocate(fzu(nd1*nx*ny,2))
allocate(fzl(nd1*nx*ny,2))
allocate(fzzu(nd2*nx*ny,2))
allocate(fzzl(nd2*nx*ny,2))
y1=0.0
y2=0.0
yf=0.0
CALL ludcmp(nx,ny,nz_all,1,0,0)
refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! Kwon
!FDM normal start==================================
if (restartnum==0) then
if (myid.eq.master) then
OPEN(305,FILE='sfield.bin',form='unformatted',status='unknown')
DO i=1,nx
READ (305) yy(i) ! Yr
ENDDO
CLOSE (305)
end if
call MPI_BCAST(yy, nx, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
in_yr=yy(1) ! inlet_Yr
out_yr=yy(nx) ! outlet_Yr
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
y1(3,ii,i,j)=0.0 ! Yr initializing
enddo
enddo
enddo
! flame location setup
sum_wrate=0.; sumc=0.
DO ii=1,nz
DO j=1,ny
DO i=1,nx
yr=y1(2,i,j,ii)/y1(1,i,j,ii)
c=1.-yr
IF (c.lt.0.) c=0.
wrate=pre*yr*exp(-ac/(1.+bc*c))
! cold boundary difficulty treatment
IF(c.le.c_ref) THEN
wrate=min_wr
IF(c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
sum_wrate=sum_wrate+wrate*(hx*hy*hy)
sumc=sumc+(1.-yr)
ENDDO
ENDDO
ENDDO
sum_wrate1=sum_wrate/(hy*hy*ny*nz_all)
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
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)))
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.
633 format (' ** Consumption Speed, Sc = ',f7.4)
634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
! Restart setup ================================================================
else
CALL fdm_restart_read
fullsavenum=restartnum+1
endif
@ -193,9 +355,12 @@ module m_fdm_calc
real*8 :: fl_location,delu,wrate,yr
real*8 :: c
real*8 :: uvel,vvel,wvel ! J. Kwon
fdmcyc=itime
! Mean velocty setup
do k=1,nz_all
do k=1,nz
do i=1,nx
do j=1,ny
u_(i,j,k)=u_(i,j,k)+dummyu_
@ -203,32 +368,34 @@ module m_fdm_calc
enddo
enddo
write(*,*)
write(*,*) '=========================================================='
if (myid.eq.0) write(*,*)
if (myid.eq.0) write(*,*) '=========================================================='
WRITE(*,932) TIME,DT
if (myid.eq.0) WRITE(*,932) TIME,DT
932 format(' Spectral results at time = ',f10.5,', dT = ',f7.5)
write(*,933) dummyu_
if (myid.eq.0) write(*,933) dummyu_
933 format(' ** Mean U = ',f7.4)
!----------------------------------------------------------------
! FDM DNS BGN
!----------------------------------------------------------------
umax=-1.0d300;umin=1.0d300;vmax=-1.0d300;vmin=1.0d300;wmax=-1.0d300;wmin=1.0d300
velmax1(1)=maxval(u_)
velmin1(1)=minval(u_)
velmax1(2)=maxval(v_)
velmin1(2)=minval(v_)
velmax1(3)=maxval(w_)
velmin1(3)=minval(w_)
!$omp parallel do private(i,j,k) reduction(max:umax,vmax,wmax) reduction(min:umin,vmin,wmin) NUM_THREADS(6) schedule(static)
do k=1,nz_all
do j=1,ny
do i=1,nx
umax=max(umax,u_(i,j,k))
umin=min(umin,u_(i,j,k))
vmax=max(vmax,v_(i,j,k))
vmin=min(vmin,v_(i,j,k))
wmax=max(wmax,w_(i,j,k))
wmin=min(wmin,w_(i,j,k))
enddo
enddo
enddo
! MPI_REDUCE [uvw](max|min)
CALL MPI_ALLREDUCE(velmax1,velmax,3,MPI_REAL8,MPI_MAX,MPI_COMM_TASK,mpi_err)
CALL MPI_ALLREDUCE(velmin1,velmin,3,MPI_REAL8,MPI_MIN,MPI_COMM_TASK,mpi_err)
umax=velmax(1)
umin=velmin(1)
vmax=velmax(2)
vmin=velmin(2)
wmax=velmax(3)
wmin=velmin(3)
uvel = max(abs(umax), abs(umin))
vvel = max(abs(vmax), abs(vmin))
@ -241,12 +408,14 @@ module m_fdm_calc
fdmdt=min(DT,visdt,convdt)
DT=fdmdt
call solve(nx,ny,nz_all,u_,v_,w_,y1,y2,yf)
call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf)
fdmtime=time+fdmdt
write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
if (myid.eq.0) write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
!$omp parallel do private(i,j,k,yr,wrate) reduction(+:sum_wrate,sumc) NUM_THREADS(6) schedule(static)
do k=1,nz_all
sumc=0.
sum_wrate=0.
do k=1,nz
do j=1,ny
do i=1,nx
yr=y1(2,i,j,k)/y1(1,i,j,k)
@ -263,21 +432,25 @@ module m_fdm_calc
enddo
enddo
! Real time results for Sc and flame location.
sum_wrate=sum_wrate/(hy*hy*REAL(ny*nz_all))
write(*,633) sum_wrate
sum_wrate1=sum_wrate/(hy*hy*REAL(ny*nz_all))
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
633 format (' ** Consumption Speed, Sc = ',f7.4)
fl_location=(hx*REAL(nx-1))*(1.-(sumc/(REAL(nx*ny*nz_all))))
sumc=sumc*(hx*hy*hy)/(hy*(REAL(ny)-1.)*hy*(REAL(ny)-1.))
write(*,634) fl_location/(REAL(nx-1)*hx)*100.
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.
634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
write(St_data,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx-1)*hx)*100.
if (myid.eq.0) write(St_data,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx-1)*hx)*100.
! Control the inflow mean velocity, dummyu_
write(*,'(a20,2f12.8)') ' ** oldsumc, sumc = ',oldsumc,sumc
if (myid.eq.0) write(*,'(a20,2f12.8)') ' ** oldsumc, sumc = ',oldsumc,sumc
if(swadtv.ne.0.and.mod((fdmcyc+1),swadtv).eq.0) then
time_int=fdmtime-time_int
delu=(sumc-oldsumc)/time_int
@ -291,11 +464,11 @@ module m_fdm_calc
! FDM DNS END
!----------------------------------------------------------------------------
write(*,*) 'Number of whole calculation = ',fdmcyc
write(*,*) '=========================================================='
write(*,*)
if (myid.eq.0) write(*,*) 'Number of whole calculation = ',fdmcyc
if (myid.eq.0) write(*,*) '=========================================================='
if (myid.eq.0) write(*,*)
call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
if (myid.eq.0) call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
return
@ -307,59 +480,172 @@ module m_fdm_calc
implicit none
integer :: i,j,k,xx,yy,zz,ii
integer :: n
integer :: idx1, idx2, idx3
real*8 :: wrate,yr,yp
real*8 :: r1_(2,xx,yy,zz),f_(2,xx,yy,zz)
real*8 :: r1_(neq,xx,yy,zz),f_(neq,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 :: ux(nd1,xx),dux(nd1,xx),d2ux(nd2,xx)
real*8 :: uy(nd1,yy),duy(nd1,yy),d2uy(nd2,yy)
real*8 :: uz(nd1,zz),duz(nd1,zz),d2uz(nd2,zz)
real*8 :: uux(nd2,xx)
real*8 :: uuy(nd2,yy)
real*8 :: uuz(nd2,zz)
real*8 :: y
! reaction source term
DO k=1,zz
DO j=1,yy
DO i=1,xx
! y=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
wrate=0.0
! wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate
! IF ((1.-y).le.c_ref) THEN
! wrate=min_wr
! 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
f_(1,i,j,k) = 0.0 ! continuity
f_(2,i,j,k) = - wrate ! species conservation
f_(3,i,j,k) = - wrate ! species conservation
ENDDO
ENDDO
ENDDO
!! 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:YF
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*YF
uz(5,k)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
uz(6,k)=uz(3,k)*uz(5,k) ! 6:rho*w*YO
uuz (1,k)=uz(2,k)
uuz (2,k)=uz(5,k)
ENDDO
CALL dfp (yy,hy, uz(1:nd1,:), duz(1:nd1,:),nd1,3)
CALL d2fp(yy,hy,uuz(1:nd2,:),d2uz(1:nd2,:),nd2,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(1,k)+duz(1,k)*duz(2,k)) ! species conserv.
f_(3,i,j,k) = f_(3,i,j,k) - duz(6,k) + diff*(uz(1,k)*d2uz(2,k)+duz(1,k)*duz(5,k)) ! species conserv.
ENDDO
ENDDO
ENDDO
else
DO k=1,zz
DO j=1,yy
DO i=1,xx
idx3 = xx*(j-1)+i
idx1 = (idx3-1)*nd1
idx2 = (idx3-1)*nd2
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:YF
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*YF
fz(idx1+5,k) = r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
fz(idx1+6,k) = r1_(3,i,j,k)*ww_(i,j,k) ! 6:rho*w*YO
fzz(idx2+1,k) = r1_(2,i,j,k)/r1_(1,i,j,k)
fzz(idx2+2,k) = r1_(3,i,j,k)/r1_(1,i,j,k)
ENDDO
ENDDO
ENDDO
CALL par_dfp (fz, fzl, fzu, dfz, hy, nd1*xx*yy, zz, yy, 3)
CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, nd2*xx*yy, zz, yy, 3)
DO k=1,zz
DO j=1,yy
DO i=1,xx
idx3 = xx*(j-1)+i
idx1 = (idx3-1)*nd1
idx2 = (idx3-1)*nd2
! -( d(rho*w)/dz )
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) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2+1,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv.
f_(3,i,j,k) = f_(3,i,j,k) - dfz(idx1+6,k) + diff*(fz(idx1+1,k)*dfzz(idx2+2,k)+dfz(idx1+1,k)*dfz(idx1+5,k)) ! species conserv.
ENDDO
ENDDO
ENDDO
endif
! x-direction
!$omp parallel do private(i,j,k,ux,dux,d2ux,wrate,uy,duy,d2uy) NUM_THREADS(6) schedule(static)
DO k=1,zz
DO j=1,yy
DO i=1,xx
ux(1,i)=r1_(1,i,j,k) ! 1:rho
ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF
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
ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*YF
ux(5,i)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
ux(6,i)=ux(3,i)*ux(5,i) ! 6:rho*u*YO
uux (1,i)=ux(2,i)
uux (2,i)=ux(5,i)
ENDDO
CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
CALL d2fnonp(xx,hx,ux(2,:),d2ux(:),1,1)
CALL dfnonp (xx,hx, ux(1:nd1,:), dux(1:nd1,:),nd1,1)
CALL d2fnonp(xx,hx,uux(1:nd2,:),d2ux(1:nd2,:),nd2,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)=-dux(3,i) ! continuity
f_(1,i,j,k) = f_(1,i,j,k) - dux(3,i) ! continuity
! -( 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)=-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(1,i)+dux(1,i)*dux(2,i)) ! species conservation
f_(3,i,j,k) = f_(3,i,j,k) - dux(6,i) + diff*(ux(1,i)*d2ux(2,i)+dux(1,i)*dux(5,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
uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF
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*YF
uy(5,j)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
uy(6,j)=uy(3,j)*uy(5,j) ! 6:rho*v*YO
uuy (1,j)=uy(2,j)
uuy (2,j)=uy(5,j)
ENDDO
CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
CALL d2fp(yy,hy,uy(2,:),d2uy(:),1,2)
CALL dfp (yy,hy, uy(1:nd1,:), duy(1:nd1,:),nd1,2)
CALL d2fp(yy,hy,uuy(1:nd2,:),d2uy(1:nd2,:),nd2,2)
DO j=1,yy
! -( d(rho*v)/dy )
@ -368,45 +654,22 @@ module m_fdm_calc
! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
! = -( d(rho*v*Yr)/dy )
! + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(j)+duy(1,j)*duy(2,j)) ! species conserv.
f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(1,j)+duy(1,j)*duy(2,j)) ! species conserv.
f_(3,i,j,k)=f_(3,i,j,k)-duy(6,j) + diff*(uy(1,j)*d2uy(2,j)+duy(1,j)*duy(5,j)) ! species conserv.
ENDDO
ENDDO
ENDDO
!! z-direction
!$omp parallel do private(i,j,k,uz,duz,d2uz) NUM_THREADS(6) schedule(static)
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
ENDDO
CALL dfp(zz,hy,uz(1:4,:),duz(1:4,:),4,3)
CALL d2fp(zz,hy,uz(2,:),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
! Boundary condition
!$omp parallel do private(k,j,i) NUM_THREADS(6) schedule(static)
DO k=1,zz
DO j=1,yy
DO i=1,yrsw
f_(2,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr
f_(3,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr
ENDDO
IF (uu_(xx,j,k).lt.0.) f_(2,xx,j,k)=f_(1,xx,j,k)*r1_(2,xx,j,k)/r1_(1,xx,j,k)
IF (uu_(xx,j,k).lt.0.) f_(3,xx,j,k)=f_(1,xx,j,k)*r1_(3,xx,j,k)/r1_(1,xx,j,k)
ENDDO
ENDDO
@ -419,7 +682,7 @@ module m_fdm_calc
integer :: istage,xx,yy,zz,i
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz)
istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
@ -437,7 +700,7 @@ module m_fdm_calc
integer :: i,j,k,xx,yy,zz
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz)
CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
@ -451,7 +714,7 @@ module m_fdm_calc
integer :: i,j,k,xx,yy,zz,istage
real*8 :: at,bt , wrate , yr
real*8 :: ri(2,xx,yy,zz),r1(2,xx,yy,zz),r2(2,xx,yy,zz),f(2,xx,yy,zz)
real*8 :: ri(neq,xx,yy,zz),r1(neq,xx,yy,zz),r2(neq,xx,yy,zz),f(neq,xx,yy,zz)
real*8 :: a(5),b(5)
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
integer :: nfinal, iscr, mspec, mpict, msave, nmindt, nv
@ -473,11 +736,10 @@ module m_fdm_calc
IF(istage<5) THEN
at=a(istage)*fdmdt
bt=(b(istage)-a(istage))*fdmdt
!$omp parallel do private(k,j,i,nv) NUM_THREADS(6) schedule(static)
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,2
DO nv=1,neq
r1(nv,i,j,k)=r1(nv,i,j,k)+at*f(nv,i,j,k)
r2(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
ENDDO
@ -486,21 +748,20 @@ module m_fdm_calc
ENDDO
ELSE
bt=b(istage)*fdmdt
sumc=0.
sum_wrate=0.
!$omp parallel do private(k,j,i,nv) NUM_THREADS(6) schedule(static)
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,2
DO nv=1,neq
r1(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
ENDDO
!==========rho=1 treatment
r1(2,i,j,k)=r1(2,i,j,k)/r1(1,i,j,k)
r1(3,i,j,k)=r1(3,i,j,k)/r1(1,i,j,k)
r1(1,i,j,k)=1.
!==========Max Yr=1 treatment
r1(2,i,j,k)=MIN(in_yr,r1(2,i,j,k))
r1(3,i,j,k)=MIN(in_yr,r1(3,i,j,k))
!==========Min Yr=0 treatment
! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))

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)

4
m_timing.f90
View file

@ -19,8 +19,8 @@ contains
! make the default job runlimit to be 6 months
job_runlimit = 6 * 30 * 24 * 60
! make the default job runlimit to be 12 hours
job_runlimit = 12 * 60
! make the default job runlimit to be 12 hours -> one month
job_runlimit = 12 * 60 * 30 * 2
write(out,'(a30,3x,i20)') 'job_runlimit (default):',job_runlimit

38
main.f90
View file

@ -67,10 +67,7 @@ program x_code
! Initializing the LES stuff
if (les) call m_les_begin
if (task.eq.'hydro' .and. myid.eq.master) then
call prepare_fdm
end if
call MPI_BARRIER(MPI_COMM_TASK, mpi_err)
call prepare_fdm
call stat_main
@ -134,20 +131,7 @@ program x_code
if (fdm_sw.ne.0) then ! Reacting
call gather_4
if (myid.eq.master) then
call fdm_exe
do id_to=1,numprocs-1
call MPI_Isend (dummy, 1, MPI_REAL, id_to, 0, MPI_COMM_TASK, mpi_request, mpi_err)
end do
else
call MPI_Irecv (dummy, 1, MPI_REAL, master, 0, MPI_COMM_TASK, mpi_request, mpi_err)
do
call MPI_Test(mpi_request, flag, mpi_status, mpi_err)
if (flag) exit
call SLEEPQQ (200)
end do
end if
call MPI_BARRIER(MPI_COMM_TASK, mpi_err)
endif
! RHS for scalars
@ -202,18 +186,18 @@ program x_code
! write the restart file if it's the time
if (mod(itime,IPRINT2).eq.0) then
call restart_write_parallel
!remove hit restart file
call system ('rm -f '//run_name//'.64.'//ex_file_ext)
if (fdm_sw.ne.0) then ! Reacting
call gather_4
if (myid.eq.master) then
!write itime, savenum pair
if(TIME.ge.bgn_save) call fdm_restart_write
write (restart_pair, *) file_ext, fullsavenum-1
call write_tp_field
!call gather_4
if (TIME.ge.bgn_save) then
call fdm_restart_write
if (myid.eq.master) then
write (restart_pair, *) file_ext, fullsavenum-1, '\n'
!call write_tp_field
end if
end if
call MPI_BARRIER(MPI_COMM_TASK, mpi_err)
endif
!remove HIT3D restart file except the latest
call system ('rm -f '//run_name//'.64.'//ex_file_ext)
ex_file_ext = file_ext
endif
@ -266,7 +250,7 @@ program x_code
if (task_split) call fields_to_stats
if (mod(itime,iprint1).eq.0) call stat_main
if (fdm_sw.eq.0) then ! non-reacting
if (mod(itime,iprint2).eq.0) call io_write_4
!if (mod(itime,iprint2).eq.0) call io_write_4
endif
end if
end if stats

4
write_stats.f90
View file

@ -31,6 +31,7 @@ subroutine write_tp_field
do i=1,nx,spx
write(tp_field,'(6e18.9)') dx*real(i),dy*real(j),dz*real(k),u_(i,j,k),v_(i,j,k),w_(i,j,k) ! u',v',w'
enddo
call flush(tp_field)
enddo
enddo
@ -49,13 +50,12 @@ subroutine write_tp_field
(1.-y1(2,i,j,ii)/y1(1,i,j,ii)),u_(i,j,ii), &
v_(i,j,ii),w_(i,j,ii)
enddo
call flush(tp_field)
enddo
enddo
endif
call flush(tp_field)
897 format('ZONE T="time= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4,' F=POINT')
end subroutine write_tp_field