research/dns-hit3d-fdm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
dns-hit3d-fdm/m_fdm_calc.f90
ignis c8d81aaf92 model spectrum written after generating velocity 
tp_field
convdt is now based on absolute maximum
sumwrate calculation is moved
omp private += i,j,k
model_e_spec
mpi_init_thread
file_units and tar_lo
kolmogorov scales
2014-05-17 15:40:19 +09:00

529 lines
16 KiB
Fortran
Raw Blame History

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
integer :: fdmcyc,fdmsavecount
real*8 :: fdmtime,fdmdt
real*8 :: sumc,oldsumc,time_int,sum_wrate
real*8 :: visdt, convdt
integer :: fdmstep
real*8 :: coe,tt1,tt2,tt3,tt4
real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon
!===========================================================================
!===========================================================================
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_,v_,w_,y1
CLOSE (fullsavenum)
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
if (fdm_sw .eq. 0) then
return
endif
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
! 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
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
!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
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
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(*,*)
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
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_all
do i=1,nx
do j=1,ny
u_(i,j,k)=u_(i,j,k)+dummyu_
enddo
enddo
enddo
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)
!----------------------------------------------------------------
! FDM DNS BGN
!----------------------------------------------------------------
umax=-1.0d300;umin=1.0d300;vmax=-1.0d300;vmin=1.0d300;wmax=-1.0d300;wmin=1.0d300
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
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_all,u_,v_,w_,y1,y2,yf)
fdmtime=time+fdmdt
write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
do k=1,nz_all
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_wrate=sum_wrate/(hy*hy*REAL(ny*nz_all))
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.
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.
! 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
time_int=fdmtime
endif
!----------------------------------------------------------------------------
! FDM DNS END
!----------------------------------------------------------------------------
write(*,*) 'Number of whole calculation = ',fdmcyc
write(*,*) '=========================================================='
write(*,*)
call write_vel_max
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(i,j,k,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,:),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)=-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,:),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(i,j,k,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,:),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(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 private(k,j,i,nv) 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.
!$omp parallel do private(k,j,i,nv) 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)+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
Reference in a new issue View git blame Copy permalink