From fd12a34f77bd68740d79ce50009577f638579886 Mon Sep 17 00:00:00 2001 From: ignis Date: Sat, 26 Apr 2014 20:16:36 +0900 Subject: [PATCH] minor changes and add new source files --- Makefile | 5 +- m_compact.f90 | 386 +++++++++++++++++++++++++++ m_fdm_calc.f90 | 680 +++++++++++++++++++++++++++++++++++++++++++++++ m_hit_result.f90 | 107 ++++++++ m_timing.f90 | 4 +- pressure.f90 | 5 +- 6 files changed, 1182 insertions(+), 5 deletions(-) create mode 100644 m_compact.f90 create mode 100644 m_fdm_calc.f90 create mode 100644 m_hit_result.f90 diff --git a/Makefile b/Makefile index 979d6c6..1a81f9d 100644 --- a/Makefile +++ b/Makefile @@ -92,7 +92,10 @@ MODULES += \ m_force.o\ m_rand_knuth.o\ m_les.o\ - RANDu.o + RANDu.o\ + m_hit_result.o\ + m_compact.o\ + m_fdm_calc.o # Objects OBJ = main.o\ diff --git a/m_compact.f90 b/m_compact.f90 new file mode 100644 index 0000000..fd175d6 --- /dev/null +++ b/m_compact.f90 @@ -0,0 +1,386 @@ + MODULE m_compact + IMPLICIT NONE + PRIVATE + REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & + lyf,lys,wyf,wys, & + lzf,lzs,wzf,wzs +! lyzf,lyzs,wyzf,wyzs + INTEGER :: nxc,nyc,nzc + REAL*8, PARAMETER :: ezero = 1.0e-14 + + PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp + + CONTAINS + + SUBROUTINE ludcmp(nxx,nyy,nzz,xp,yp,zp) + INTEGER, INTENT(IN) :: nxx,nyy,nzz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nxx + nyc=nyy + nzc=nzz +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + ALLOCATE(lxf(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lxs(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(xp.eq.0) THEN + ALLOCATE(wxf(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wxs(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + CALL p_lud(1,nxc) + ELSE + CALL nonp_lud(1,nxc) + ENDIF + + ALLOCATE(lyf(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lys(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(yp.eq.0) THEN + ALLOCATE(wyf(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wys(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + CALL p_lud(2,nyc) + ELSE + call nonp_lud(2,nyc) + ENDIF + + ALLOCATE(lzf(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lzs(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(zp.eq.0) THEN + ALLOCATE(wzf(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wzs(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + CALL p_lud(3,nzc) + ELSE + call nonp_lud(3,nzc) + ENDIF + +! CALL x_lud +! CALL yz_lud + + END SUBROUTINE ludcmp + + SUBROUTINE nonp_lud(xyz,xx) + INTEGER :: i,xyz,xx + REAL*8, DIMENSION(xx) :: aa + aa=3. + aa(1)=0.5 ; aa(2)=4. + aa(xx-1)=4. ; aa(xx)=0.5 +! first derivative + IF (xyz.eq.1) CALL stdlu(aa,xx,lxf) ! x-direction + IF (xyz.eq.2) CALL stdlu(aa,xx,lyf) ! y-direction + IF (xyz.eq.3) CALL stdlu(aa,xx,lzf) ! z-direction + aa=5.5 + aa(1)=2./11. ; aa(2)=10. + aa(xx-1)=10. ; aa(xx)=2./11. +! second derivative + IF (xyz.eq.1) CALL stdlu(aa,xx,lxs) ! x-direction + IF (xyz.eq.2) CALL stdlu(aa,xx,lys) ! y-direction + IF (xyz.eq.3) CALL stdlu(aa,xx,lzs) ! z-direction + END SUBROUTINE nonp_lud + + + SUBROUTINE p_lud(xyz,xx) + INTEGER :: i,xyz,xx + REAL*8 :: a + a=3. ! first derivative + IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction + IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction + IF (xyz.eq.3) CALL ptdlu(a,xx,lzf,wzf) ! z-direction + a=11./2. ! second derivative + IF (xyz.eq.1) CALL ptdlu(a,xx,lxs,wxs) ! x-direction + IF (xyz.eq.2) CALL ptdlu(a,xx,lys,wys) ! y-direction + IF (xyz.eq.3) CALL ptdlu(a,xx,lzs,wzs) ! z-direction + END SUBROUTINE p_lud + + SUBROUTINE stdlu(a,n,l) + INTEGER :: n + REAL*8 :: a(n),l(n) + REAL*8 :: d + INTEGER :: i + l(1)=1.0/a(1) + DO i=2,n + d=a(i)-l(i-1) + l(i)=1.0/d + ENDDO + END SUBROUTINE stdlu + + SUBROUTINE ptdlu(a,n,l,w) + INTEGER :: n + REAL*8 :: a,l(n),w(n) + INTEGER :: i + REAL*8 :: aa(n),d + + DO i=1,n-1 + aa(i)=a + ENDDO + i=n-1 + call stdlu(aa,i,l) + w(1)=1.0 + DO i=2,n-2 + w(i)=-l(i-1)*w(i-1) + ENDDO + w(n-1)=1.0-l(n-2)*w(n-2) + DO i=1,n-1 + w(i)=w(i)*l(i) + ENDDO + d=a + DO i=1,n-1 + d=d-w(i)*w(i)/l(i) + ENDDO + l(n)=1./d + END SUBROUTINE ptdlu + + SUBROUTINE dfnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: 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,r3,a,b,c,h1,t1,t2,t3,t4 + + h1=1./h + + r1=7./3. + r2=1./12. + r3=3. + a=-1.25 + b=1. + c=0.25 + + DO j=1,nd + dx(j,n-1)=x(j,n)-x(j,n-2) + dx(j,n)=-(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)) + dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)) + dx(j,2)=x(j,3)-x(j,1) + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n)=0. + IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3)) dx(j,1)=0. + dx(j,n-1)=dx(j,n-1)*h1*r3 + dx(j,n)=dx(j,n)*h1 + dx(j,1)=dx(j,1)*h1 + dx(j,2)=dx(j,2)*h1*r3 + ENDDO + + DO i=3,n-2 + DO j=1,nd + t1=x(j,i+1)-x(j,i-1) + t2=x(j,i+2)-x(j,i-2) + dx(j,i)=h1*(r1*t1+r2*t2) + ENDDO + ENDDO + IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction + IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction + IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction + END SUBROUTINE dfnonp + + SUBROUTINE dfp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: 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 + + h1=1./h + r1=7./3. + r2=1./12. + + DO j=1,nd + dx(j,n-1)=(r1*(x(j,n)-x(j,n-2))+r2*(x(j,1)-x(j,n-3))) + dx(j,n)=(r1*(x(j,1)-x(j,n-1))+r2*(x(j,2)-x(j,n-2))) + dx(j,1)=(r1*(x(j,2)-x(j,n))+r2*(x(j,3)-x(j,n-1))) + dx(j,2)=(r1*(x(j,3)-x(j,1))+r2*(x(j,4)-x(j,n))) + dx(j,n-1)=dx(j,n-1)*h1 + dx(j,n)=dx(j,n)*h1 + dx(j,1)=dx(j,1)*h1 + dx(j,2)=dx(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))) + dx(j,i)=dx(j,i)*h1 + ENDDO + ENDDO + + IF (dir.eq.1) CALL ptdslv(dx,n,lxf,wxf,nd) ! x-direction + IF (dir.eq.2) CALL ptdslv(dx,n,lyf,wyf,nd) ! y-direction + IF (dir.eq.3) CALL ptdslv(dx,n,lzf,wzf,nd) ! z-direction + + END SUBROUTINE dfp + + SUBROUTINE ptdslv(r,n,l,w,nd) + INTEGER,INTENT(IN) :: n,nd + REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL*8,INTENT(IN),DIMENSION(:) :: l,w + INTEGER i,j + REAL*8, DIMENSION(nd) :: sum + DO j=1,nd + sum(j)=w(1)*r(j,1) + r(j,1)=r(j,1)*l(1) + ENDDO + DO i=2,n-1 + DO j=1,nd + 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=1,nd + r(j,n)=l(n)*(r(j,n)-sum(j)) + r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n) + ENDDO + DO i=n-2,1,-1 + DO j=1,nd + r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n) + ENDDO + ENDDO + END SUBROUTINE ptdslv + + SUBROUTINE d2fp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: 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 + h2=1./(h*h) + r1=6. + r2=3./8. + DO j=1,nd + t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2)) + t2 = (x(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 (x(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) + + t1 = (x(j,1)-2.*x(j,n)+x(j,n-1)) + t2 = (x(j,2)-2.*x(j,n)+x(j,n-2)) + IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0. + IF (x(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) + + t1 = (x(j,2)-2.*x(j,1)+x(j,n)) + t2 = (x(j,3)-2.*x(j,1)+x(j,n-1)) + IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0. + IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-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) + + t1 = (x(j,3)-2.*x(j,2)+x(j,1)) + t2 = (x(j,4)-2.*x(j,2)+x(j,n)) + 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.x(j,n)) 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) + + dx(j,n-1)=dx(j,n-1)*h2 + dx(j,n)=dx(j,n)*h2 + dx(j,1)=dx(j,1)*h2 + dx(j,2)=dx(j,2)*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) + dx(j,i)=dx(j,i)*h2 + ENDDO + ENDDO + IF (dir.eq.1) CALL ptdslv(dx,n,lxs,wxs,nd) ! x-direction + IF (dir.eq.2) CALL ptdslv(dx,n,lys,wys,nd) ! y-direction + IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction + END SUBROUTINE d2fp + + SUBROUTINE tdslv(r,n,l,nd) + INTEGER,INTENT(IN) :: n,nd + REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL*8,INTENT(IN),DIMENSION(:) :: l + INTEGER i,j + REAL*8 t1 + DO j=1,nd + r(j,1)=r(j,1)*l(1) + ENDDO + DO i=2,n + DO j=1,nd + t1=r(j,i)-r(j,i-1) + r(j,i)=l(i)*t1 + ENDDO + ENDDO + DO i=n-1,1,-1 + DO j=1,nd + r(j,i)=r(j,i)-l(i)*r(j,i+1) + ENDDO + ENDDO + END SUBROUTINE tdslv + + SUBROUTINE d2fnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: 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,r3,a,b,c,e,t1,t2 + + h2=1./(h*h) + r1=6. + r2=3./8. + r3=12. + a=13./11. + b=-27./11. + c=15./11. + e=-1./11. + + DO j=1,nd + dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)+e*x(j,4)) + dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1)) + dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2)) + dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3)) + IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0. + IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0. + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0. + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0. + dx(j,1)=dx(j,1)*h2 + dx(j,2)=dx(j,2)*h2*r3 + dx(j,n-1)=dx(j,n-1)*h2*r3 + dx(j,n)=dx(j,n)*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) + dx(j,i)=dx(j,i)*h2 + ENDDO + ENDDO + + IF (dir.eq.1) CALL tdslv(dx,n,lxs,nd) ! x-direction + IF (dir.eq.2) CALL tdslv(dx,n,lys,nd) ! y-direction + IF (dir.eq.3) CALL tdslv(dx,n,lzs,nd) ! z-direction + + END SUBROUTINE d2fnonp + + END MODULE m_compact diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 new file mode 100644 index 0000000..33218bb --- /dev/null +++ b/m_fdm_calc.f90 @@ -0,0 +1,680 @@ +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. ! + avgr=0. ! + + 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) ! + avgr=avgr/REAL(ny*nz) ! + + 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) ! + avgr=avgr/REAL(ny*nz) ! + + 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 + + diff --git a/m_hit_result.f90 b/m_hit_result.f90 new file mode 100644 index 0000000..4795dd9 --- /dev/null +++ b/m_hit_result.f90 @@ -0,0 +1,107 @@ +!================================================================================ +! m_hit_result - module for printing out the results +! such as energy spectrum and velocity field +! +! Time-stamp: <2014-03-19 J. Kwon> +!================================================================================ + +module m_hit_result + +contains + + subroutine result_files_open + use m_parameters +! statistics + write(*,*) 'result_open,ITMIN',ITMIN + + if(ITMIN.eq.0) then + open(69,file='tp_stat.dat',form='formatted') + write(69,*) 'VARIABLES = "itime","time","tke","rms_u`","int_LS","Eddy_TO_time","Kol_LS"' ! + write(69,*) '"Kol_VS","Kol_TS","Taylor_LS","dissipation"' + else + open(69,file='tp_stat.dat',status='old',position='append') + endif + +! flucating velocity field, u',v',w' + write(*,*) 'result_files_open, fdm_sw',fdm_sw + if(fdm_sw.eq.0) then + if(ITMIN.eq.0) then + open(101,FILE='tp_field.dat') + write(101,*) 'VARIABLES = "x","y","z","u","v","w"' + else + open(101,file='tp_field.dat',status='old',position='append') + endif + else + if(ITMIN.eq.0) then + open(101,file='tp_field.dat') + write(101,*) 'VARIABLES = "X","Y","Z","C","U","V","W"' + else + open(101,file='tp_field.dat',status='old',position='append') + endif + endif + +! model spectrum + if(ITMIN.eq.0) then + open(unit=404,file="model_spectrum.dat") + write(404,*) 'VARIABLES = "k","Model_E(k)"' ! + write(404,*) '"Kolmogorov","Exponential","VonKarman"' ! + endif + +! realtime energy/dissipation spectrum + if(ITMIN.eq.0) then + open(unit=501,file="tp_spec.dat") + write(501,*) 'VARIABLES = "k","E(k)","D(k)"' + else + open(unit=501,file="tp_spec.dat",status='old',position='append') + endif + +! realtime energy spectrum + if(ITMIN.eq.0) then + open(unit=502,file="tp_tke_spec.dat") + write(502,*) 'VARIABLES = "k","E(k)"' + else + open(unit=502,file="tp_tke_spec.dat",status='old',position='append') + endif +! realtime dissipation spectrum + if(ITMIN.eq.0) then + open(unit=503,file="tp_eps_spec.dat") + write(503,*) 'VARIABLES = "k","D(k)"' + else + open(unit=503,file="tp_eps_spec.dat",status='old',position='append') + endif + +! realtime max, min velocities + if(ITMIN.eq.0) then + open(unit=504,file="tp_max_min_uvw.dat") + write(504,*) 'VARIABLES = "time","max_u","max_v","max_w","min_u","min_v","min_w"' + else + open(unit=504,file="tp_max_min_uvw.dat",status='old',position='append') + endif +! FDM temporary output +! 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"," + if(fdm_sw.ne.0) then + if(ITMIN.eq.0) then + open(505,file='St_data.dat') + write(505,*)'VARIABLES = "Time","Mean Velocity","Sc","Flame Location"' + else + open(505,file='St_data.dat',status='old',position='append') + endif + endif + +! if(fdm_sw.ne.0) then +! if(ITMIN.eq.0) then +! open(506,file='3D_field.dat') +! write(506,*) 'VARIABLES = "X","Y","Z","C","U","V","W"' +! else +! open(506,file='3D_field.dat',status='old',position='append') +! endif +! endif + + + + end subroutine result_files_open + +end module m_hit_result diff --git a/m_timing.f90 b/m_timing.f90 index a627a54..d043117 100644 --- a/m_timing.f90 +++ b/m_timing.f90 @@ -22,7 +22,7 @@ contains ! make the default job runlimit to be 12 hours job_runlimit = 12 * 60 - write(out,*) 'job_runlimit (default):',job_runlimit + write(out,'(a30,3x,i20)') 'job_runlimit (default):',job_runlimit if(myid_world.eq.0) then @@ -35,7 +35,7 @@ contains end if call MPI_BCAST(job_runlimit,1,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err) - write(out,*) 'job_runlimit: ',job_runlimit + write(out,'(a30,3x,i20)') 'job_runlimit:',job_runlimit call flush(out) return diff --git a/pressure.f90 b/pressure.f90 index 4298cb6..4944968 100644 --- a/pressure.f90 +++ b/pressure.f90 @@ -84,8 +84,9 @@ subroutine divergence call MPI_REDUCE(d1,dmax,1,MPI_REAL8,MPI_MAX,0,MPI_COMM_TASK,mpi_err) if (myid.eq.0) then - write(out,*) 'divergence:',dmin,dmax -!!$ print *, 'divergence:',dmin,dmax + write(out,*) + write(out,'(a30,3x,2e20.10)') 'divergence: min, max',dmin,dmax + write(out,*) call flush(out) end if