720 lines
24 KiB
Fortran
720 lines
24 KiB
Fortran
module m_fdm_calc
|
|
|
|
use m_parameters
|
|
use m_io
|
|
use m_compact
|
|
|
|
implicit none
|
|
|
|
!variables
|
|
|
|
real*8, dimension(:,:,:), allocatable :: u_,v_,w_
|
|
real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf
|
|
|
|
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
|
|
|
|
integer :: vtype1, vtype2
|
|
integer, dimension(128) :: scnt1, sdisp1, stype1, rcnt1, rdisp1, rtype1
|
|
integer, dimension(128) :: scnt2, sdisp2, stype2, rcnt2, rdisp2, rtype2
|
|
|
|
|
|
!===========================================================================
|
|
!===========================================================================
|
|
|
|
contains
|
|
|
|
|
|
subroutine fdm_restart_read
|
|
integer(kind=MPI_INTEGER_KIND) :: fh
|
|
integer(kind=MPI_OFFSET_KIND) :: offset
|
|
|
|
write(fname,"('fort.',i4)") restartnum
|
|
|
|
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,file=fname,form='unformatted',status='unknown')
|
|
read (restartnum) fdmtime,tt1,tt2,tt3,oldsumc,time_int
|
|
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
|
|
CLOSE (restartnum)
|
|
|
|
if(fdmtime.ge.t_fdmsave) t_fdmsave=t_fdmsave+dt_fdmsave
|
|
if(fdmtime.ge.t_fullsave) t_fullsave=t_fullsave+dt_fullsave
|
|
|
|
write (6,456) fdmtime,fdmcyc
|
|
456 format(' Restart time : ',f10.5,' / FDM cycle : ',i6)
|
|
write(*,*)
|
|
write(*,*) ' Save Times : FDM FULL '
|
|
write(*,454) t_fdmsave,t_fullsave
|
|
454 format(' ',f7.3,' ',f7.3)
|
|
write(*,*)' Save Intervals : FDM FULL','from input'
|
|
write(*,453) dt_fdmsave,dt_fullsave
|
|
453 format(' ',f7.3,' ',f7.3)
|
|
write(*,*)
|
|
|
|
if (ABS(fdmtime-TIME).le.1.0e-10) then
|
|
write(*,*) ' Spectral and FDM times are consistent.'
|
|
else
|
|
write(*,*)' !! Warning : Spectral and FDM times are different !!'
|
|
write(*,455) TIME,fdmtime
|
|
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)
|
|
|
|
count = nx * ny * nz
|
|
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 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 = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 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 = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 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 = 2 * nx * ny * nz
|
|
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 3*nx*ny*nz_all*8 + myid*2*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
|
|
! --------------- 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
|
|
|
|
|
|
count = nx * ny * nz
|
|
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 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 = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 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 = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 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 = 2 * nx * ny * nz
|
|
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 3*nx*ny*nz_all*8 + myid*2*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(2,nx,ny,nz))
|
|
allocate(y2(2,nx,ny,nz))
|
|
allocate(yf(2,nx,ny,nz))
|
|
|
|
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
|
|
|
|
|
|
|
|
!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 i=1,ny
|
|
do j=1,nz
|
|
do ii=1,nx
|
|
y1(1,ii,i,j)=1. ! rho initializing
|
|
y1(2,ii,i,j)=yy(ii) ! 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
|
|
|
|
fl_location=(hx*(nx-1.))*(1.-(sumc/(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
|
|
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
|
|
|
|
! End of Restart setup ========================================================
|
|
|
|
return
|
|
end subroutine prepare_fdm
|
|
|
|
subroutine fdm_exe
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k,ii
|
|
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
|
|
do i=1,nx
|
|
do j=1,ny
|
|
u_(i,j,k)=u_(i,j,k)+dummyu_
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
if (myid.eq.0) write(*,*)
|
|
if (myid.eq.0) write(*,*) '=========================================================='
|
|
|
|
if (myid.eq.0) WRITE(*,932) TIME,DT
|
|
932 format(' Spectral results at time = ',f10.5,', dT = ',f7.5)
|
|
if (myid.eq.0) write(*,933) dummyu_
|
|
933 format(' ** Mean U = ',f7.4)
|
|
|
|
!----------------------------------------------------------------
|
|
! FDM DNS BGN
|
|
!----------------------------------------------------------------
|
|
velmax1(1)=maxval(u_)
|
|
velmin1(1)=minval(u_)
|
|
velmax1(2)=maxval(v_)
|
|
velmin1(2)=minval(v_)
|
|
velmax1(3)=maxval(w_)
|
|
velmin1(3)=minval(w_)
|
|
|
|
! 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))
|
|
wvel = max(abs(wmax), abs(wmin))
|
|
umax = umax - dummyu_
|
|
umin = umin - dummyu_
|
|
|
|
visdt=max(0.,0.3*fdmcfl*hx**2./nu)
|
|
convdt=max(0.,fdmcfl*hx/(uvel+vvel+wvel))
|
|
|
|
fdmdt=min(DT,visdt,convdt)
|
|
DT=fdmdt
|
|
call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf)
|
|
fdmtime=time+fdmdt
|
|
if (myid.eq.0) write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
|
|
|
|
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)
|
|
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))
|
|
IF((1.-yr).le.c_ref) THEN
|
|
wrate=min_wr
|
|
IF((1.-yr).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-yr-c_ref))+ &
|
|
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
|
ENDIF
|
|
!get sum_wrate
|
|
sum_wrate=sum_wrate+wrate*(hx*hy*hy)
|
|
sumc=sumc+(1.-yr)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
|
|
! Real time results for Sc and flame location.
|
|
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))))
|
|
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.')
|
|
|
|
|
|
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_
|
|
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
|
|
if ((fdmtime).gt.startad.and.oldsumc.ne.0.) then
|
|
dummyu_=dummyu_+delu
|
|
endif
|
|
time_int=fdmtime
|
|
endif
|
|
|
|
!----------------------------------------------------------------------------
|
|
! FDM DNS END
|
|
!----------------------------------------------------------------------------
|
|
|
|
if (myid.eq.0) write(*,*) 'Number of whole calculation = ',fdmcyc
|
|
if (myid.eq.0) write(*,*) '=========================================================='
|
|
if (myid.eq.0) write(*,*)
|
|
|
|
if (myid.eq.0) call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
|
|
|
|
return
|
|
|
|
end subroutine fdm_exe
|
|
|
|
|
|
SUBROUTINE fns(r1_,f_,xx,yy,zz,uu_,vv_,ww_)
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k,xx,yy,zz,ii
|
|
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)
|
|
|
|
|
|
!CALL MPI_ALLTOALLW (ww_, scnt1, sdisp1, stype1, ww_t, rcnt1, rdisp1, rtype1, MPI_COMM_TASK, mpi_err)
|
|
|
|
!CALL MPI_ALLTOALLW (r1_, scnt2, sdisp2, stype2, r1_t, rcnt2, rdisp2, rtype2, MPI_COMM_TASK, mpi_err)
|
|
|
|
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
|
|
|
|
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
|
|
|
|
!! 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
|
|
|
|
! x-direction
|
|
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(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
|
|
ENDDO
|
|
|
|
CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
|
|
CALL d2fnonp(xx,hx,ux(2,:),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
|
|
|
|
! -( 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
|
|
|
|
ENDDO
|
|
ENDDO
|
|
|
|
!! y-direction
|
|
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
|
|
ENDDO
|
|
|
|
CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
|
|
CALL d2fp(yy,hy,uy(2,:),d2uy(:),1,2)
|
|
|
|
DO j=1,yy
|
|
! -( d(rho*v)/dy )
|
|
f_(1,i,j,k)=f_(1,i,j,k)-duy(3,j) ! continuity
|
|
|
|
! -( 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.
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
|
|
|
|
! Boundary condition
|
|
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
|
|
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)
|
|
ENDDO
|
|
ENDDO
|
|
|
|
return
|
|
END SUBROUTINE fns
|
|
|
|
subroutine RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
|
|
|
implicit none
|
|
|
|
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)
|
|
|
|
|
|
istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
|
|
istage=2; CALL substep(yy1,yy2,yy1,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
|
|
istage=3; CALL substep(yy2,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
|
|
istage=4; CALL substep(yy1,yy2,yy1,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
|
|
istage=5; CALL substep(yy2,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
|
|
|
|
return
|
|
END SUBROUTINE RK4
|
|
|
|
SUBROUTINE solve(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
|
|
|
IMPLICIT NONE
|
|
|
|
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)
|
|
|
|
|
|
CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
|
|
|
|
return
|
|
END SUBROUTINE solve
|
|
|
|
subroutine substep(ri,r1,r2,f,xx,yy,zz,istage,uu_,vv_,ww_)
|
|
|
|
implicit none
|
|
|
|
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 :: 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
|
|
|
|
|
|
a(1)= 970286171893./4311952581923.
|
|
a(2)= 6584761158862./12103376702013.
|
|
a(3)= 2251764453980./15575788980749.
|
|
a(4)= 26877169314380./34165994151039.
|
|
a(5)=0.
|
|
b(1)= 1153189308089./22510343858157.
|
|
b(2)= 1772645290293./4653164025191.
|
|
b(3)= -1672844663538./4480602732383.
|
|
b(4)= 2114624349019./3568978502595.
|
|
b(5)= 5198255086312./14908931495163.
|
|
|
|
CALL fns(ri,f,xx,yy,zz,uu_,vv_,ww_)
|
|
|
|
IF(istage<5) THEN
|
|
at=a(istage)*fdmdt
|
|
bt=(b(istage)-a(istage))*fdmdt
|
|
DO k=1,zz
|
|
DO j=1,yy
|
|
DO i=1,xx
|
|
DO nv=1,2
|
|
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
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ELSE
|
|
bt=b(istage)*fdmdt
|
|
DO k=1,zz
|
|
DO j=1,yy
|
|
DO i=1,xx
|
|
DO nv=1,2
|
|
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(1,i,j,k)=1.
|
|
!==========Max Yr=1 treatment
|
|
r1(2,i,j,k)=MIN(in_yr,r1(2,i,j,k))
|
|
!==========Min Yr=0 treatment
|
|
! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))
|
|
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
|
|
return
|
|
END SUBROUTINE substep
|
|
|
|
|
|
end module m_fdm_calc
|
|
|
|
|