diff --git a/Makefile b/Makefile index f1005a9..f68f052 100644 --- a/Makefile +++ b/Makefile @@ -115,7 +115,8 @@ OBJ = main.o\ velocity_rescale.o\ write_tmp4.o\ gather_tmp4.o\ - gather_4.o + gather_4.o\ + write_stats.o # ------------------------------------------------------- diff --git a/gather_tmp4.f90 b/gather_tmp4.f90 index b6aaf0f..1263352 100644 --- a/gather_tmp4.f90 +++ b/gather_tmp4.f90 @@ -16,13 +16,13 @@ subroutine gather_tmp4(u_i) if (myid.ne.master) then id_to = master tag = myid - call MPI_SEND(tmp4,count,MPI_REAL4,master,tag,MPI_COMM_TASK,mpi_err) + call MPI_SEND(tmp4,count,MPI_REAL8,master,tag,MPI_COMM_TASK,mpi_err) else u_i(:,:,1:nz) = tmp4(:,:,:) do id_from=1,numprocs-1 tag = id_from - call MPI_RECV(tmp4,count,MPI_REAL4,id_from,tag,MPI_COMM_TASK,mpi_status,mpi_err) + call MPI_RECV(tmp4,count,MPI_REAL8,id_from,tag,MPI_COMM_TASK,mpi_status,mpi_err) u_i(:,:,id_from*nz+1:id_from*nz+nz) = tmp4(:,:,:) end do end if diff --git a/init_velocity.f90 b/init_velocity.f90 index f9e4f5f..1a5cbf7 100644 --- a/init_velocity.f90 +++ b/init_velocity.f90 @@ -209,13 +209,6 @@ subroutine init_velocity ! by L. Guichard, J. Reveillon and R. Hauguel ! ------------------------------------------------------------------------ -! open(unit=404,file="model_spectrum.dat") -! write(404,*) 'VARIABLES = "X","Energy(k)"' ! -! write(404,*) '"Kolmogorov","Exponential","VonKarman"' ! - -! write(502,'(a15,f11.5,a1)')'ZONE T= " time=',time,'"' -! write(502,*) 'I=',kmax,'F=POINT' - if (isp_type.eq.4) then !ksj model_e_spec=zip k0=(SQRT(two*PI))/t_length @@ -223,72 +216,8 @@ subroutine init_velocity do i=1,kmax k=real(i,8) model_e_spec(i)=model_A*((k/k0)**4.0)*exp(-2.0*((k/k0)**2.0)) -! e_spec_log00=log10(model_e_spec(i)) - - wmag=real(k, 8) - ratio=wmag / peak_wavenum - e_spec_log01=wmag**(-5.d0/3.d0) - e_spec_log02=ratio**3 / peak_wavenum * exp(-3.0D0*ratio) - - fac = two * PI * ratio - e_spec_log03=fac**4 / (one + fac**2)**3 - - if(model_e_spec(i).le.0.) e_spec_log00=0. -! write(404,'(i3,4e15.5)') i,model_e_spec(i),& -! e_spec_log01,e_spec_log02,e_spec_log03 -! write(502,'(i5,e15.5)') i,model_e_spec(i) enddo -! close(404) - -! write target stats at t=0 on tp_stats.dat - ! getting the total energy - model_tke = sum(model_e_spec(1:kmax)) - - ! finding dissipation spectrum and total dissipation - do k = 1,kmax - model_d_spec(k) = model_e_spec(k) * real(k**2,8) * two * nu - end do - model_eps_v = sum(model_d_spec(1:kmax)) - - ! finding Kolmogorov scale - model_eta = (nu**3/model_eps_v)**0.25 ! Kolmogorov length scale - - ! Kolmogorov time and velocity scale, J. Kwon - kol_ts=sqrt(nu/model_eps_v) ! Kolmogorov time scale(kol_ts) - kol_vs=model_eta/kol_ts ! Kolmogorov velocity scale(kol_vs) - - model_uvar=two/three*model_tke - ! rms velocity fluctuation, J. Kwon - rms_u_prime=sqrt(model_uvar) - - ! integral length scale(x_length) - model_sctmp = zip - do k = 1, kmax - model_sctmp = model_sctmp + model_e_spec(k) / real(k,8) - end do - model_x_length = PI / two * model_sctmp / model_uvar - - ! Taylor microscale - model_lambda = sqrt(15.d0 * model_uvar * nu / model_eps_v) - - ! Eddy turnover time - model_tau_e = model_x_length / sqrt(model_uvar) - - - ! outputting all this in the stat1 file -! inquire(file='tp_stat.dat', exist=there, opened=there2) -! if (.not.there) 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"' -! end if -! if(there.and..not.there2) then -! open(69,file='tp_stat.dat',position='append') -! end if -! write(69,"(i10,f15.10,9e15.6)") itime,time,model_tke,rms_u_prime,model_x_length, & -! model_tau_e,model_eta,kol_vs,kol_ts,model_lambda,& -! model_eps_v endif @@ -350,6 +279,7 @@ subroutine init_velocity write(out,*) "Generated the velocities." call flush(out) + call write_model_spectrum ! deallocate work arrays deallocate(e_spec, e_spec1, rr, hits, hits1, stat=ierr) diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 5c0e175..25e4467 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -18,6 +18,7 @@ module m_fdm_calc real*8 :: visdt, convdt integer :: fdmstep real*8 :: coe,tt1,tt2,tt3,tt4 + real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon !=========================================================================== !=========================================================================== @@ -46,6 +47,7 @@ module m_fdm_calc integer :: i,j,ii,k real*8 :: fl_location,wrate,yr,c + real*8 :: tar_sumc if (fdm_sw .eq. 0) then return @@ -78,6 +80,7 @@ module m_fdm_calc 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 @@ -143,6 +146,7 @@ module m_fdm_calc 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) @@ -203,7 +207,7 @@ module m_fdm_calc integer :: i,j,k,ii real*8 :: fl_location,delu,wrate,yr real*8 :: c - real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon + real*8 :: uvel,vvel,wvel ! J. Kwon fdmcyc=itime ! Mean velocty setup do k=1,nz_all @@ -239,30 +243,38 @@ module m_fdm_calc 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/(umax+vmax+wmax)) - - WRITE(101,897) fdmtime,svfx,svfy,svfy - - do ii=1,nz_all !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 + convdt=max(0.,fdmcfl*hx/(uvel+vvel+wvel)) fdmdt=min(DT,visdt,convdt) - call solve(nx,ny,nz_all,u_,v_,w_,y1,y2,yf) 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 @@ -274,7 +286,7 @@ module m_fdm_calc 634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.') - write(505,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx)*hx)*100. + 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 @@ -283,35 +295,9 @@ module m_fdm_calc if ((fdmtime).gt.startad.and.oldsumc.ne.0.) then dummyu_=dummyu_-delu endif - oldsumc=sumc time_int=fdmtime endif -! 3D_field.dat ================================================================ - tp_field:if(mod(itime,iprint1).eq.0) then - write(6,*) ' ## 3D Field is being written.' - WRITE(101,897) fdmtime,svfx,svfy,svfy -897 format('ZONE T="time= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4) - - do ii=1,nz_all !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 - - fdmsavecount=fdmsavecount+1 - t_fdmsave=t_fdmsave+dt_fdmsave - endif tp_field -!================================================ - - write(504,'(f12.6,6e18.8)')time,umax,vmax,wmax,umin,vmin,wmin - !---------------------------------------------------------------------------- ! FDM DNS END !---------------------------------------------------------------------------- @@ -320,6 +306,8 @@ module m_fdm_calc write(*,*) '==========================================================' write(*,*) + call write_vel_max + return end subroutine fdm_exe @@ -339,7 +327,7 @@ module m_fdm_calc real*8 :: Ly(yy,zz),Dy ! x-direction -!$omp parallel do private(ux,dux,d2ux,wrate,uy,duy,d2uy) NUM_THREADS(4) schedule(static) +!$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 @@ -397,7 +385,7 @@ module m_fdm_calc ENDDO !! z-direction -!$omp parallel do private(uz,duz,d2uz) NUM_THREADS(4) schedule(static) +!$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 @@ -423,7 +411,7 @@ module m_fdm_calc ENDDO ! Boundary condition -!$omp parallel do NUM_THREADS(4) schedule(static) +!$omp parallel do private(k,j,i) NUM_THREADS(4) schedule(static) DO k=1,zz DO j=1,yy DO i=1,yrsw @@ -496,7 +484,7 @@ module m_fdm_calc IF(istage<5) THEN at=a(istage)*fdmdt bt=(b(istage)-a(istage))*fdmdt -!$omp parallel do NUM_THREADS(4) schedule(static) +!$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 @@ -511,29 +499,22 @@ module m_fdm_calc 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)) -!========= - yr=r1(2,i,j,k)/r1(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) +!==========Min Yr=0 treatment +! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k)) + ENDDO ENDDO ENDDO diff --git a/m_force.f90 b/m_force.f90 index f14e993..db69d0d 100644 --- a/m_force.f90 +++ b/m_force.f90 @@ -86,7 +86,6 @@ contains end do ! Restarting setup - if(isp_type.eq.4) then allocate(model_e_spec(kmax)) ; model_e_spec=zip k0=(SQRT(two*PI))/t_length model_A=16.0*(SQRT(two/PI))*(t_u_prime**2.0)/k0 @@ -94,7 +93,6 @@ contains k=real(i,8) model_e_spec(i)=model_A*((k/k0)**4.0)*exp(-2.0*((k/k0)**2.0)) enddo - endif !!$ ! writing out the nodes !!$ do n = 1,n_forced_nodes diff --git a/m_hit_result.f90 b/m_hit_result.f90 index 4795dd9..9a84588 100644 --- a/m_hit_result.f90 +++ b/m_hit_result.f90 @@ -6,6 +6,8 @@ !================================================================================ module m_hit_result + use m_parameters + implicit none contains @@ -15,67 +17,67 @@ contains 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"' + open(tp_stat,file='tp_stat.dat',form='formatted') + write(tp_stat,*) 'VARIABLES = "itime","time","tke","rms_u`","int_LS","Eddy_TO_time","Kol_LS"' ! + write(tp_stat,*) '"Kol_VS","Kol_TS","Taylor_LS","dissipation"' else - open(69,file='tp_stat.dat',status='old',position='append') + open(tp_stat,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"' + open(tp_field,FILE='tp_field.dat') + write(tp_field,*) 'VARIABLES = "x","y","z","u","v","w"' else - open(101,file='tp_field.dat',status='old',position='append') + open(tp_field,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"' + open(tp_field,file='tp_field.dat') + write(tp_field,*) 'VARIABLES = "X","Y","Z","C","U","V","W"' else - open(101,file='tp_field.dat',status='old',position='append') + open(tp_field,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"' ! + open(unit=model_spectrum,file="model_spectrum.dat") + write(model_spectrum,*) 'VARIABLES = "k","Model_E(k)"' ! + write(model_spectrum,*) '"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)"' + open(unit=tp_spec,file="tp_spec.dat") + write(tp_spec,*) 'VARIABLES = "k","E(k)","D(k)"' else - open(unit=501,file="tp_spec.dat",status='old',position='append') + open(unit=tp_spec,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)"' + open(unit=tp_tke_spec,file="tp_tke_spec.dat") + write(tp_tke_spec,*) 'VARIABLES = "k","E(k)"' else - open(unit=502,file="tp_tke_spec.dat",status='old',position='append') + open(unit=tp_tke_spec,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)"' + open(unit=tp_eps_spec,file="tp_eps_spec.dat") + write(tp_eps_spec,*) 'VARIABLES = "k","D(k)"' else - open(unit=503,file="tp_eps_spec.dat",status='old',position='append') + open(unit=tp_eps_spec,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"' + open(unit=tp_max_min_uvw,file="tp_max_min_uvw.dat") + write(tp_max_min_uvw,*) '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') + open(unit=tp_max_min_uvw,file="tp_max_min_uvw.dat",status='old',position='append') endif ! FDM temporary output ! OPEN(23,FILE='St_data.dat') @@ -84,10 +86,10 @@ contains ! 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"' + open(St_data,file='St_data.dat') + write(St_data,*)'VARIABLES = "Time","Mean Velocity","Sc","Flame Location"' else - open(505,file='St_data.dat',status='old',position='append') + open(St_data,file='St_data.dat',status='old',position='append') endif endif diff --git a/m_mpi.f90 b/m_mpi.f90 index d0c8fee..3cc221a 100644 --- a/m_mpi.f90 +++ b/m_mpi.f90 @@ -20,7 +20,7 @@ module m_openmpi logical :: iammaster integer(kind=MPI_INTEGER_KIND) :: myid_world, numprocs_world integer(kind=MPI_INTEGER_KIND) :: numprocs_hydro, numprocs_stats, numprocs_parts - integer(kind=MPI_INTEGER_KIND) :: myid, numprocs, master, mpi_err, mpi_info + integer(kind=MPI_INTEGER_KIND) :: myid, numprocs, master, mpi_err, mpi_info, mpi_provide integer(kind=MPI_INTEGER_KIND) :: id_to, id_from, tag, count integer(kind=MPI_INTEGER_KIND) :: id_root_hydro, id_root_stats, id_root_parts @@ -64,7 +64,7 @@ contains ! initializing MPI environment - call MPI_INIT(mpi_err) + call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err) call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err) call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err) diff --git a/m_openmpi.f90 b/m_openmpi.f90 index 4ae1f51..73d297a 100644 --- a/m_openmpi.f90 +++ b/m_openmpi.f90 @@ -20,7 +20,7 @@ module m_openmpi logical :: iammaster integer(kind=MPI_INTEGER_KIND) :: myid_world, numprocs_world integer(kind=MPI_INTEGER_KIND) :: numprocs_hydro, numprocs_stats, numprocs_parts - integer(kind=MPI_INTEGER_KIND) :: myid, numprocs, master, mpi_err, mpi_info + integer(kind=MPI_INTEGER_KIND) :: myid, numprocs, master, mpi_err, mpi_info, mpi_provide integer(kind=MPI_INTEGER_KIND) :: id_to, id_from, tag, count integer(kind=MPI_INTEGER_KIND) :: id_root_hydro, id_root_stats, id_root_parts @@ -64,7 +64,7 @@ contains ! initializing MPI environment - call MPI_INIT(mpi_err) + call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err) call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err) call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err) diff --git a/m_parameters.f90 b/m_parameters.f90 index 0aadf64..6c26934 100644 --- a/m_parameters.f90 +++ b/m_parameters.f90 @@ -114,6 +114,19 @@ module m_parameters real*8 :: ac,bc,c_cut,c_ref,min_wr,prof_wr,hx,hy,diff,pre real*8 :: dt_fullsave,dt_fdmsave,t_fullsave,t_fdmsave real*8 :: fdmcfl,startad,dummyu_,schmidt + real*8 :: tar_lo + +! files for statistics + integer, parameter :: tp_stat = 690 + integer, parameter :: tp_field = 101 + integer, parameter :: tp_plane = 506 + integer, parameter :: model_spectrum = 404 + integer, parameter :: tp_spec = 501 + integer, parameter :: tp_tke_spec = 502 + integer, parameter :: tp_eps_spec = 503 + integer, parameter :: tp_max_min_uvw = 504 + integer, parameter :: St_data = 505 + !================================================================================ contains @@ -566,6 +579,8 @@ contains read(in,*,ERR=9000,END=9000) fdmcfl write(out,'(a30,3x,e20.10)') 'FDM_CFL =',fdmcfl + read(in,*,ERR=9000,END=9000) tar_lo + write(out,'(a30,3x,e20.10)') 'tar_lo =',tar_lo ! initialize ---------------------- t_fdmsave=dt_fdmsave diff --git a/m_stats.f90 b/m_stats.f90 index b651cb2..e7c5fc8 100644 --- a/m_stats.f90 +++ b/m_stats.f90 @@ -8,6 +8,7 @@ module m_stats real*8 :: energy, eps_v, eta, etakmax, enstrophy, re_lambda, uvar, x_length real*8 :: lambda, re_lambda1, tau_e real*8 :: sctmp + real*8 :: kol_ts,kol_vs,rms_u_prime contains @@ -125,12 +126,20 @@ contains eps_v = sum(e_spec1(1:kmax)) ! finding Kolmogorov scale - eta = (nu**3/eps_v)**0.25 - etakmax = eta * real(kmax,8) + eta = (nu**3/eps_v)**0.25 ! Kolmogorov length scale + etakmax = eta * real(kmax,8) ! Kolmogorov length scale * max_wavenum + + ! Kolmogorov time and velocity scale, J. Kwon + kol_ts=sqrt(nu/eps_v) ! Kolmogorov time scale(kol_ts) + kol_vs=eta/kol_ts ! Kolmogorov velocity scale(kol_vs) ! variance uvar = two/three*energy - ! integral length scale + + ! rms velocity fluctuation, J. Kwon + rms_u_prime=sqrt(uvar) + + ! integral length scale(x_length) sctmp = zip do k = 1, kmax sctmp = sctmp + e_spec(k) / real(k,8) @@ -152,7 +161,7 @@ contains inquire(file='stat1.gp', exist=there, opened=there2) if (.not.there) then open(69,file='stat1.gp',form='formatted') - write(69,*) '# 1.itime 2.time 3.energy 4.diss 5.eta 6.enstrophy 7.R_lambda' + write(69,'(A)') '# 1.itime 2.time 3.energy 4.diss 5.eta 6.enstrophy 7.R_lambda' end if if(there.and..not.there2) then open(69,file='stat1.gp',position='append') @@ -165,7 +174,7 @@ contains inquire(file='stat2.gp', exist=there, opened=there2) if (.not.there) then open(70,file='stat2.gp',form='formatted') - write(70,'(A)') '# 1.itime 2.time 3.int LS 4. lambda 5.R_lambda1 6.tau_e 7.etakmax' + write(70,'(A)') '# 1.itime 2.time 3.int_LS 4.lambda 5.R_lambda1 6.tau_e 7.etakmax' end if if(there.and..not.there2) then open(70,file='stat2.gp',position='append') @@ -184,6 +193,8 @@ contains end do close(900) + call write_tp_stat + end if return end subroutine stat_velocity diff --git a/main.f90 b/main.f90 index 2ead4d6..7303ec1 100644 --- a/main.f90 +++ b/main.f90 @@ -51,6 +51,7 @@ program x_code ! getting the wallclock runlimit for the job call get_job_runlimit + call result_files_open !----------------------------------------------------------------------- ! Starting from the beginning or from the saved flowfield @@ -64,13 +65,12 @@ program x_code ! Initializing the LES stuff if (les) call m_les_begin - !call result_files_open - if (fdm_sw.ne.0) then if (task.eq.'hydro' .and. myid.eq.master) then call prepare_fdm end if call MPI_BARRIER(MPI_COMM_TASK, mpi_err) - end if + + call stat_main ! checking divergence if (task.eq.'hydro') call divergence diff --git a/write_stats.f90 b/write_stats.f90 new file mode 100644 index 0000000..56b651f --- /dev/null +++ b/write_stats.f90 @@ -0,0 +1,110 @@ +subroutine write_tp_stat + use m_parameters + use m_stats + implicit none + + write(tp_stat,"(i10,f15.10,9e15.6)") itime,time, & + energy,rms_u_prime,x_length,tau_e,eta,kol_vs, & + kol_ts,lambda,eps_v + call flush(tp_stat) + +end subroutine write_tp_stat + + +subroutine write_tp_field + use m_parameters + use m_fdm_calc + implicit none + + integer :: i, j, k, ii + + if (fdm_sw.eq.0) then + write(tp_field,897) time,nx/spx,ny/spy,nz/spz + do k=1,nz_all,spz + do j=1,ny,spy + do i=1,nx,spx + write(tp_field,'(6e18.9)') dx*real(i),dy*real(j),dz*real(k),u_(i,j,k),v_(i,j,k),w_(i,j,k) ! u',v',w' + enddo + enddo + enddo + else + write(tp_field,897) time,nx/svf,ny/svf,nz_all/svf + + do ii=1,nz_all,svf + do j=1,ny,svf + do i=1,nx,svf + WRITE(tp_field,'(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) + enddo + enddo + enddo + + endif + + call flush(tp_field) + +897 format('ZONE T="time= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4,' F=POINT') + +end subroutine write_tp_field + +subroutine write_model_spectrum + use m_parameters + implicit none + integer :: i + real*8 :: k, wmag, ratio, fac + real*8 :: e_spec_log00=0. + real*8 :: e_spec_log01=0. + real*8 :: e_spec_log02=0. + real*8 :: e_spec_log03=0. + + do i=1,kmax + k=real(i,8) + + wmag=real(k, 8) + ratio=wmag / peak_wavenum + e_spec_log01=wmag**(-5.d0/3.d0) + e_spec_log02=ratio**3 / peak_wavenum * exp(-3.0D0*ratio) + + fac = two * PI * ratio + e_spec_log03=fac**4 / (one + fac**2)**3 + + if(model_e_spec(i).le.0.) e_spec_log00=0. + write(model_spectrum,'(i3,4e15.5)') i,model_e_spec(i),& + e_spec_log01,e_spec_log02,e_spec_log03 + write(tp_tke_spec,'(i5,e15.5)') i,model_e_spec(i) + enddo +end subroutine write_model_spectrum + +subroutine write_vel_max + use m_parameters + use m_fdm_calc + implicit none + write(tp_max_min_uvw,'(f12.6,6e18.8)')time,umax,vmax,wmax,umin,vmin,wmin +end subroutine write_vel_max + +subroutine write_tp_eps_spec + use m_parameters + use m_stats + implicit none + integer :: i + write(tp_eps_spec,'(a15,f11.5,a1)')'ZONE T= " time=',time,'"' + write(tp_eps_spec,*) 'I=',kmax,'F=POINT' + do i=1,kmax + write(tp_eps_spec,'(i5,e15.5)') i,e_spec1(i) + enddo +end subroutine write_tp_eps_spec + +subroutine write_tp_spec + use m_parameters + use m_stats + implicit none + integer :: i + + write(501,'(a15,f11.5,a1)')'ZONE T= " time=',time,'"' + write(501,*) 'I=',kmax,'F=POINT' + do i=1,kmax + write(501,'(i5,2e15.5)') i,e_spec(i),e_spec1(i) + enddo +end subroutine write_tp_spec + diff --git a/write_tmp4.f90 b/write_tmp4.f90 index d0e5f7c..6525d17 100644 --- a/write_tmp4.f90 +++ b/write_tmp4.f90 @@ -17,7 +17,7 @@ subroutine write_tmp4 if (myid.ne.master) then id_to = master tag = myid - call MPI_SEND(tmp4,count,MPI_REAL4,master,tag,MPI_COMM_TASK,mpi_err) + call MPI_SEND(tmp4,count,MPI_REAL8,master,tag,MPI_COMM_TASK,mpi_err) else !! open(my_out,file=fname,form='binary') open(my_out,file=fname,form='unformatted', access='stream') @@ -26,7 +26,7 @@ subroutine write_tmp4 do id_from=1,numprocs-1 tag = id_from - call MPI_RECV(tmp4,count,MPI_REAL4,id_from,tag,MPI_COMM_TASK,mpi_status,mpi_err) + call MPI_RECV(tmp4,count,MPI_REAL8,id_from,tag,MPI_COMM_TASK,mpi_status,mpi_err) write(my_out) (((tmp4(i,j,k),i=1,nx),j=1,ny),k=1,nz) end do close(my_out)