680 lines
21 KiB
Fortran
680 lines
21 KiB
Fortran
module m_fdm_calc
|
|
|
|
use m_parameters
|
|
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 !,rest_sw
|
|
integer :: fdmcyc,fdmsavecount
|
|
real*8 :: fdmtime,fdmdt
|
|
real*8 :: sumc,oldsumc,time_int,sum_wrate
|
|
real*8 :: max_u,max_v,max_w,min_u,min_v,min_w
|
|
real*8 :: visdt, convdt
|
|
integer :: fdmstep
|
|
|
|
!===========================================================================
|
|
!===========================================================================
|
|
|
|
contains
|
|
|
|
|
|
subroutine prepare_fdm
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,ii,k
|
|
real*8 :: avgc(nx),avgr(nx)
|
|
|
|
allocate(u_(nx,ny,nz))
|
|
allocate(v_(nx,ny,nz))
|
|
allocate(w_(nx,ny,nz))
|
|
|
|
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
|
|
! rest_sw=0
|
|
|
|
allocate(y1(2,nx,ny,nz))
|
|
allocate(y2(2,nx,ny,nz))
|
|
allocate(yf(2,nx,ny,nz))
|
|
|
|
CALL ludcmp(nx,ny,nz,1,0,0)
|
|
|
|
!FDM normal start==================================
|
|
|
|
! OPEN(44,FILE='max_min_u.dat')
|
|
! write(44,*)'VARIABLES = "Time","max_u","max_v","max_w","min_u","min_v","min_w"'
|
|
! OPEN(23,FILE='St_data.dat')
|
|
! write(23,*)'VARIABLES = "Time","Mean Velocity","Sc","Flame Location"'
|
|
! OPEN(22,FILE='3D_Field.dat')
|
|
! WRITE(22,*) 'VARIABLES = "X","Y","Z","C","U","V","W"'
|
|
|
|
! refwr=pre*1.*exp(-ac/(1.+bc*ccut))
|
|
refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! Kwon
|
|
! minf=exp((c_ref-ccut)*prof_wr)
|
|
if(svf.eq.0) then
|
|
spx=1
|
|
spy=1
|
|
svfx=nx !nx
|
|
svfy=ny !ny
|
|
else
|
|
svfx=0
|
|
svfy=0
|
|
do i=2,nx,spx
|
|
svfx=svfx+1 ! # of points in x-dir. in 3D_field.dat
|
|
enddo
|
|
do i=2,ny,spy
|
|
svfy=svfy+1 ! # of points in y-dir. in 3D_field.dat
|
|
enddo
|
|
endif
|
|
|
|
|
|
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
|
|
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
|
|
|
|
|
|
! WRITE(22,*) 'VARIABLES = "X","Y","Z","C","U","V","W"'
|
|
! WRITE(22,897) fdmtime,svfx,svfy,svfy
|
|
! WRITE(101,897) fdmtime,svfx,svfy,svfy
|
|
!897 format('ZONE T="TIME= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4)
|
|
|
|
|
|
avgc=0. ! <c>
|
|
avgr=0. ! <rho>
|
|
|
|
do k=1,nz !k
|
|
do j=1,ny
|
|
do i=1,nx
|
|
|
|
avgc(i)=avgc(i)+(1.-y1(2,i,j,k)/y1(1,i,j,k))
|
|
avgr(i)=avgr(i)+y1(1,i,j,k)
|
|
|
|
! if (mod(k,spy).eq.svf.and.mod(j,spy).eq.svf.and.mod(i,spx).eq.svf) then
|
|
! WRITE(22,'(7e30.20)')REAL(i)*hx,REAL(j)*hx,REAL(k)*hx,&
|
|
! (1.-y1(2,i,j,k)/y1(1,i,j,k)),u_(i,j,k),v_(i,j,k),w_(i,j,k)
|
|
! WRITE(101,'(7e30.20)')REAL(i)*hx,REAL(j)*hx,REAL(k)*hx,&
|
|
! (1.-y1(2,i,j,k)/y1(1,i,j,k)),u_(i,j,k),v_(i,j,k),w_(i,j,k)
|
|
! endif
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
avgc=avgc/REAL(ny*nz) ! <c>
|
|
avgr=avgr/REAL(ny*nz) ! <rho>
|
|
|
|
endif
|
|
|
|
!FDM normal start======================================
|
|
|
|
return
|
|
end subroutine prepare_fdm
|
|
|
|
subroutine fdm_exe
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k,ii
|
|
real*8 :: coe,tt1,tt2,tt3,tt4
|
|
real*8 :: avgc(nx),avgr(nx),fl_location,delu,wrate,yr
|
|
real*8 :: c
|
|
real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon
|
|
! Mean velocty setup
|
|
do k=1,nz
|
|
do i=1,nx
|
|
do j=1,ny
|
|
u_(i,j,k)=u_(i,j,k)+dummyu_
|
|
! IF (i.le.inx1) THEN
|
|
! u_(i,j,k)=dummyu_
|
|
! v_(i,j,k)=0.
|
|
! w_(i,j,k)=0.
|
|
! ELSE IF (i.le.(inx1+inx2)) THEN
|
|
! coe=sin(0.5/REAL(inx2)*REAL(i-inx1-1)*PI)
|
|
! u_(i,j,k)=dummyu_+u_(i,j,k)*coe
|
|
! v_(i,j,k)=v_(i,j,k)*coe
|
|
! w_(i,j,k)=w_(i,j,k)*coe
|
|
! ENDIF
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
! Restart setup ================================================================
|
|
if (restartnum.ne.0) then
|
|
458 format(a55)
|
|
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')
|
|
read (restartnum) fdmtime,tt1,tt2,tt3,oldsumc,time_int
|
|
read (restartnum) fdmcyc,tt4,dummyu_
|
|
read (restartnum) tt1,tt2 !dt_fdmsave,dt_fullsave
|
|
read (restartnum) t_fdmsave,t_fullsave
|
|
read (restartnum) in_yr,out_yr
|
|
fdmdt=DT
|
|
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
|
|
|
|
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(*,*)
|
|
|
|
TIME=TIME-DT
|
|
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
|
|
TIME=TIME+DT
|
|
|
|
! rest_sw=0
|
|
|
|
fullsavenum=restartnum+1
|
|
restartnum=0
|
|
avgc=0.; avgr=0.
|
|
do ii=1,nz
|
|
do j=1,ny
|
|
do i=1,nx
|
|
avgc(i)=avgc(i)+(1.-y1(2,i,j,ii)/y1(1,i,j,ii))
|
|
avgr(i)=avgr(i)+y1(1,i,j,ii)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
avgc=avgc/REAL(ny*nz) ! <c>
|
|
avgr=avgr/REAL(ny*nz) ! <rho>
|
|
|
|
endif
|
|
! End of Restart setup ========================================================
|
|
|
|
write(*,*)
|
|
write(*,*) '=========================================================='
|
|
|
|
WRITE(*,932) TIME,DT
|
|
932 format(' Spectral results at time = ',f10.5,', dT = ',f7.5)
|
|
write(*,933) dummyu_
|
|
933 format(' ** Mean U = ',f7.4)
|
|
|
|
! if(itime.eq.1.or.rest_sw.eq.1) then
|
|
|
|
! flame location setup
|
|
if(itime.eq.1) then
|
|
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))
|
|
|
|
! IF ((1.-yr).le.ccut) THEN
|
|
! wrate=pre*yr*exp(-c_ref/(1.+bc*(1.-yr)))
|
|
! ENDIF
|
|
|
|
! 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)
|
|
write(*,633) sum_wrate
|
|
fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz)))
|
|
sumc=sumc*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.))
|
|
write(*,634) fl_location/(REAL(nx-1)*hx)*100.
|
|
|
|
endif
|
|
|
|
!----------------------------------------------------------------
|
|
! FDM DNS BGN
|
|
!----------------------------------------------------------------
|
|
! fdm:if (itime.ne.1) then
|
|
! if (fdmcyc.eq.1) then
|
|
if (fdmcyc.eq.0) then
|
|
umax=0.; umin=0.; vmax=0.; vmin=0.; wmax=0.; wmin=0.
|
|
do k=1,nz
|
|
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
|
|
! convdt=DT/2.
|
|
! visdt=convdt
|
|
visdt=max(0.,0.3*fdmcfl*hx**2./nu)
|
|
convdt=max(0.,fdmcfl*hx/(umax+vmax+wmax))
|
|
|
|
WRITE(101,897) fdmtime,svfx,svfy,svfy
|
|
!897 format('ZONE T="TIME= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4)
|
|
do ii=1,nz !k
|
|
do j=1,ny
|
|
do i=1,nx
|
|
if (mod(ii,spy).eq.svf.and.mod(j,spy).eq.svf.and.mod(i,spx).eq.svf) then
|
|
WRITE(101,'(7e30.20)') REAL(i)*hx,REAL(j)*hx,REAL(ii)*hx, &
|
|
(1.-y1(2,i,j,ii)/y1(1,i,j,ii)),u_(i,j,ii), &
|
|
v_(i,j,ii),w_(i,j,ii)
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
endif
|
|
|
|
! DQ's fdmdt setup
|
|
! do ii=10,1,-1
|
|
! tt1=DT/REAL(ii)
|
|
! if (tt1.le.MIN(convdt,visdt)) then
|
|
! fdmdt=tt1
|
|
! fdmstep=ii
|
|
! endif
|
|
! enddo
|
|
|
|
! do ii=1,fdmstep ! Time marching loop
|
|
! IF (ii.eq.fdmstep) fdmdt=TIME-fdmtime
|
|
|
|
! fdmtime=fdmtime+fdmdt
|
|
! write(*,931) ii,fdmstep,fdmtime,fdmdt,MIN(convdt,visdt)
|
|
! call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf)
|
|
!931 format(' FDM : ',i2,'/',i2,', time = ',f10.5,', FDM dT = ',f7.5,' < ',f7.5)
|
|
! enddo
|
|
! DQ's fdmdt setup
|
|
|
|
fdmdt=min(DT,visdt,convdt)
|
|
call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf)
|
|
DT=fdmdt
|
|
fdmtime=fdmtime+fdmdt
|
|
write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
|
|
! Real time results for Sc and flame location.
|
|
sum_wrate=sum_wrate/(hy*hy*REAL(ny*nz))
|
|
write(*,633) sum_wrate
|
|
633 format (' ** Consumption Speed, Sc = ',f7.4)
|
|
fl_location=(hx*REAL(nx-1))*(1.-(sumc/(REAL(nx*ny*nz))))
|
|
sumc=sumc*(hx*hy*hy)/(hy*(REAL(ny)-1.)*hy*(REAL(ny)-1.))
|
|
write(*,634) fl_location/(REAL(nx-1)*hx)*100.
|
|
634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
|
|
|
|
! write(23,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx)*hx)*100.
|
|
write(505,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx)*hx)*100.
|
|
|
|
! Control the inflow mean velocity, dummyu_
|
|
if(swadtv.ne.0.and.mod((fdmcyc+1),swadtv).eq.0) then
|
|
time_int=fdmtime-time_int
|
|
delu=(oldsumc-sumc)/time_int
|
|
if ((fdmtime).gt.startad.and.oldsumc.ne.0.) then
|
|
dummyu_=dummyu_-delu
|
|
endif
|
|
oldsumc=sumc
|
|
time_int=fdmtime
|
|
endif
|
|
|
|
! Full save ===================================================================
|
|
if(fdmtime.ge.t_fullsave) then
|
|
write(*,*) '======================================================='
|
|
write(*,*) 'Full results are being written',fullsavenum
|
|
OPEN (fullsavenum,form='unformatted',status='unknown')
|
|
write (fullsavenum) fdmtime,nx,ny,nz,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_,v_,w_,y1
|
|
CLOSE (fullsavenum)
|
|
write(*,*) '======================================================='
|
|
fullsavenum=fullsavenum+1
|
|
t_fullsave=t_fullsave+dt_fullsave
|
|
endif
|
|
|
|
! 3D_field.dat ================================================================
|
|
tp_field:if(fdmtime.ge.t_fdmsave) then
|
|
write(6,*) ' ## 3D Field is being written.'
|
|
! OPEN(22,FILE='3D_Field.dat')
|
|
! WRITE(22,*) 'VARIABLES = "X","Y","Z","C","U","V","W"'
|
|
! WRITE(22,897) fdmtime,svfx,svfy,svfy
|
|
WRITE(101,897) fdmtime,svfx,svfy,svfy
|
|
897 format('ZONE T="time= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4)
|
|
|
|
do ii=1,nz !k
|
|
do j=1,ny
|
|
do i=1,nx
|
|
if (mod(ii,spy).eq.svf.and.mod(j,spy).eq.svf.and.mod(i,spx).eq.svf) then
|
|
! WRITE(22,'(7e30.20)') REAL(i)*hx,REAL(j)*hx,REAL(ii)*hx, &
|
|
! (1.-y1(2,i,j,ii)/y1(1,i,j,ii)),u_(i,j,ii), &
|
|
! v_(i,j,ii),w_(i,j,ii)
|
|
WRITE(101,'(7e30.20)') REAL(i)*hx,REAL(j)*hx,REAL(ii)*hx, &
|
|
(1.-y1(2,i,j,ii)/y1(1,i,j,ii)),u_(i,j,ii), &
|
|
v_(i,j,ii),w_(i,j,ii)
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
fdmsavecount=fdmsavecount+1
|
|
t_fdmsave=t_fdmsave+dt_fdmsave
|
|
endif tp_field
|
|
!================================================
|
|
|
|
visdt=max(0.,0.3*fdmcfl*hx**2./nu)
|
|
convdt=max(0.,fdmcfl*hx/(max_u+max_v+max_w))
|
|
|
|
! write(44,'(e14.4,6e18.8)') TIME,max_u,max_v,max_w,min_u,min_v,min_w
|
|
if(fdm_sw.ne.0) write(504,'(f12.6,6e18.8)')time,max_u,max_v,max_w,min_u,min_v,min_w
|
|
|
|
! endif fdm
|
|
!----------------------------------------------------------------------------
|
|
! FDM DNS END
|
|
!----------------------------------------------------------------------------
|
|
|
|
write(*,*) 'Number of whole calculation = ',fdmcyc
|
|
write(*,*) '=========================================================='
|
|
fdmcyc=fdmcyc+1
|
|
write(*,*)
|
|
|
|
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
|
|
|
|
! x-direction
|
|
!$omp parallel do private(ux,dux,d2ux,wrate,uy,duy,d2uy) NUM_THREADS(4) 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(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:2,:),d2ux(:),1,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.ccut) THEN
|
|
! wrate=pre*ux(2,i)*exp(-c_ref/(1.+bc*(1.-ux(2,i)))) !wrate
|
|
! ENDIF
|
|
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
|
|
|
|
! -( 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
|
|
|
|
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:2,:),d2uy(:),1,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
|
|
|
|
!! z-direction
|
|
!$omp parallel do private(uz,duz,d2uz) NUM_THREADS(4) 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:2,:),d2uz(:),1,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 NUM_THREADS(4) 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
|
|
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
|
|
!$omp parallel do NUM_THREADS(4) schedule(static)
|
|
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
|
|
sumc=0.
|
|
sum_wrate=0.
|
|
max_u=0.; max_v=0.; max_w=0.
|
|
min_u=0.; min_v=0.; min_w=0.
|
|
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))
|
|
!=========
|
|
yr=r1(2,i,j,k)/r1(1,i,j,k)
|
|
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))
|
|
! IF ((1.-yr).le.ccut) THEN
|
|
! wrate=pre*yr*exp(-c_ref/(1.+bc*(1.-yr)))
|
|
! ENDIF
|
|
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
|
|
sum_wrate=sum_wrate+wrate*(hx*hy*hy)
|
|
sumc=sumc+(1.-yr)
|
|
!get sum_wrate
|
|
max_u=MAX(max_u,uu_(i,j,k))
|
|
max_v=MAX(max_v,vv_(i,j,k))
|
|
max_w=MAX(max_w,ww_(i,j,k))
|
|
min_u=MIN(min_u,uu_(i,j,k))
|
|
min_v=MIN(min_v,vv_(i,j,k))
|
|
min_w=MIN(min_w,ww_(i,j,k))
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
|
|
return
|
|
END SUBROUTINE substep
|
|
|
|
|
|
end module m_fdm_calc
|
|
|
|
|