diff --git a/Makefile b/Makefile index 1a81f9d..f1005a9 100644 --- a/Makefile +++ b/Makefile @@ -113,7 +113,9 @@ OBJ = main.o\ rhs_velocity.o\ rhs_scalars.o\ velocity_rescale.o\ - write_tmp4.o + write_tmp4.o\ + gather_tmp4.o\ + gather_4.o # ------------------------------------------------------- @@ -133,5 +135,5 @@ $(OBJ): $(MODULES) $(MPIF90) $(FCFLAGS) $< clean: - rm *.o *.mod $(PROG) + rm -f *.o *.mod $(PROG) diff --git a/gather_4.f90 b/gather_4.f90 new file mode 100644 index 0000000..b53fecd --- /dev/null +++ b/gather_4.f90 @@ -0,0 +1,93 @@ +subroutine gather_4 + + ! Writing out the velocities and scalars in X-space + ! to the real*4 file + + use m_parameters + use m_fields + use m_work + use x_fftw + use m_les + use m_fdm_calc + + implicit none + + integer :: n_out, n, i, j, k + real*8 :: wmag2, rkmax2 + + ! every variable will undergo a mode truncation for all modes + ! that are higher than kmax. This will ensure that the written + ! variables are isotropic + + rkmax2 = real(kmax,8)**2 + + ! number of variables to write out + n_out = 3 + if (int_scalars) n_out = n_out + n_scalars + if (les .and. n_les>0) n_out = n_out + n_les + + ! putting all variables in wrk array + do k = 1,nz + do j = 1,ny + do i = 1,nx+2 + wmag2 = akx(i)**2 + aky(k)**2 + akz(j)**2 + + if (wmag2 .gt. rkmax2) then + wrk(i,j,k,1:n_out) = zip + else + wrk(i,j,k,1:n_out) = fields(i,j,k,1:n_out) + end if + + end do + end do + end do + + ! velocities + call xFFT3d(-1,1) + fname = 'u.'//file_ext + tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1) + call gather_tmp4(u_) + + call xFFT3d(-1,2) + fname = 'v.'//file_ext + tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2) + call gather_tmp4(v_) + + call xFFT3d(-1,3) + fname = 'w.'//file_ext + tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3) + call gather_tmp4(w_) + + ! scalars + if (int_scalars) then + do n = 1,n_scalars + call xFFT3d(-1,3+n) + write(fname,"('sc',i2.2,'.',a6)") n,file_ext + tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3+n) + call write_tmp4 + + end do + end if + + ! LES quantities + if (les) then + ! turbulent viscosity + if (allocated(turb_visc)) then + write(fname,"('nu_t.',a6)") file_ext + tmp4 = turb_visc + call write_tmp4 + end if + + if (n_les > 0) then + do n = 1, n_les + call xFFT3d(-1,3+n_scalars+n) + write(fname,"('les',i1,'.',a6)") n,file_ext + tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3+n_scalars+n) + call write_tmp4 + end do + end if + + end if + + return +end subroutine gather_4 diff --git a/gather_tmp4.f90 b/gather_tmp4.f90 new file mode 100644 index 0000000..b6aaf0f --- /dev/null +++ b/gather_tmp4.f90 @@ -0,0 +1,33 @@ +subroutine gather_tmp4(u_i) + + use m_parameters + use m_io + use m_work + implicit none + + integer :: i,j,k + real*8, dimension(nx,ny,nz_all) :: u_i + +!====================================================================== + + ! --- defining the size of whole array + count = nx*ny*nz + + if (myid.ne.master) then + id_to = master + tag = myid + call MPI_SEND(tmp4,count,MPI_REAL4,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) + u_i(:,:,id_from*nz+1:id_from*nz+nz) = tmp4(:,:,:) + end do + end if + +!====================================================================== + + return +end subroutine gather_tmp4 diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 33218bb..1f44145 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -32,9 +32,9 @@ module m_fdm_calc 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)) + allocate(u_(nx,ny,nz_all)) + allocate(v_(nx,ny,nz_all)) + allocate(w_(nx,ny,nz_all)) u_=0.0 v_=0.0 @@ -57,11 +57,11 @@ module m_fdm_calc 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)) + 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,1,0,0) + CALL ludcmp(nx,ny,nz_all,1,0,0) !FDM normal start================================== @@ -104,7 +104,7 @@ module m_fdm_calc out_yr=y1(2,nx,1,1) ! outlet_Yr do i=1,ny - do j=1,nz + 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 @@ -122,7 +122,7 @@ module m_fdm_calc avgc=0. ! avgr=0. ! - do k=1,nz !k + do k=1,nz_all !k do j=1,ny do i=1,nx @@ -140,8 +140,8 @@ module m_fdm_calc enddo enddo - avgc=avgc/REAL(ny*nz) ! - avgr=avgr/REAL(ny*nz) ! + avgc=avgc/REAL(ny*nz_all) ! + avgr=avgr/REAL(ny*nz_all) ! endif @@ -160,7 +160,7 @@ module m_fdm_calc real*8 :: c real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon ! Mean velocty setup - do k=1,nz + do k=1,nz_all do i=1,nx do j=1,ny u_(i,j,k)=u_(i,j,k)+dummyu_ @@ -225,7 +225,7 @@ module m_fdm_calc fullsavenum=restartnum+1 restartnum=0 avgc=0.; avgr=0. - do ii=1,nz + do ii=1,nz_all do j=1,ny do i=1,nx avgc(i)=avgc(i)+(1.-y1(2,i,j,ii)/y1(1,i,j,ii)) @@ -234,8 +234,8 @@ module m_fdm_calc enddo enddo - avgc=avgc/REAL(ny*nz) ! - avgr=avgr/REAL(ny*nz) ! + avgc=avgc/REAL(ny*nz_all) ! + avgr=avgr/REAL(ny*nz_all) ! endif ! End of Restart setup ======================================================== @@ -253,7 +253,7 @@ module m_fdm_calc ! flame location setup if(itime.eq.1) then sum_wrate=0.; sumc=0. - DO ii=1,nz + DO ii=1,nz_all DO j=1,ny DO i=1,nx yr=y1(2,i,j,ii)/y1(1,i,j,ii) @@ -278,9 +278,9 @@ module m_fdm_calc ENDDO ENDDO - sum_wrate=sum_wrate/(hy*hy*ny*nz) + sum_wrate=sum_wrate/(hy*hy*ny*nz_all) write(*,633) sum_wrate - fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz))) + fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz_all))) sumc=sumc*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.)) write(*,634) fl_location/(REAL(nx-1)*hx)*100. @@ -293,7 +293,7 @@ module m_fdm_calc ! 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 k=1,nz_all do j=1,ny do i=1,nx umax=max(umax,u_(i,j,k)) @@ -312,7 +312,7 @@ module m_fdm_calc 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 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 @@ -339,21 +339,21 @@ module m_fdm_calc ! fdmtime=fdmtime+fdmdt ! write(*,931) ii,fdmstep,fdmtime,fdmdt,MIN(convdt,visdt) -! call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf) +! call solve(nx,ny,nz_all,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) + call solve(nx,ny,nz_all,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)) + 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)))) + 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.') @@ -377,7 +377,7 @@ module m_fdm_calc 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) 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 @@ -398,7 +398,7 @@ module m_fdm_calc 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 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 diff --git a/m_work.f90 b/m_work.f90 index c6f34cf..79d02cf 100644 --- a/m_work.f90 +++ b/m_work.f90 @@ -13,7 +13,7 @@ module m_work implicit none ! --- work arrays - real*4,allocatable :: tmp4(:,:,:) + real*8, allocatable :: tmp4(:,:,:) real*8, allocatable :: wrk(:,:,:,:) diff --git a/main.f90 b/main.f90 index 68f7632..7033025 100644 --- a/main.f90 +++ b/main.f90 @@ -13,6 +13,8 @@ program x_code use m_particles use m_filter_xfftw use m_les + use m_hit_result ! J. Kwon + use m_fdm_calc implicit none @@ -34,6 +36,10 @@ program x_code call m_stats_init call m_force_init + !call result_files_open + if (task.eq.'hydro' .and. myid.eq.master) then + call prepare_fdm + end if ! allocating and initializing particles if (task.eq.'parts') then @@ -220,8 +226,15 @@ program x_code ! these are executed regardless of the processor configuration if (task_split) call fields_to_stats if (mod(itime,iprint1).eq.0) call stat_main - if (mod(itime,iwrite4).eq.0) call io_write_4 - + if (fdm_sw.eq.0) then ! non-reacting + if (mod(itime,iwrite4).eq.0) call io_write_4 + else ! reacting + call gather_4 + if (task.eq.'hydro' .and. myid.eq.master) then + call fdm_exe + end if + call MPI_BARRIER(MPI_COMM_TASK, mpi_err) + endif end if end if stats