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_+dummyu_,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 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 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 return 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 !$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 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 !$omp parallel do private(i,j,k,yr,wrate) reduction(+:sum_wrate,sumc) NUM_THREADS(6) schedule(static) 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_ 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 !---------------------------------------------------------------------------- write(*,*) 'Number of whole calculation = ',fdmcyc write(*,*) '==========================================================' write(*,*) 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 ! 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(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(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 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(6) 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(6) 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