research/dns-hit3d-fdm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

velocity field (u,v,w) is gathered and passed to fdm part

Browse source
This commit is contained in:
ignis 2014-04-29 20:28:58 +09:00
parent abdc5f5fdc
commit db40cf48b4
6 changed files with 172 additions and 31 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
Makefile
View file

@ -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)

93
gather_4.f90 Normal file
View file

@ -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

33
gather_tmp4.f90 Normal file
View file

@ -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

52
m_fdm_calc.f90
View file

@ -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. ! <c>
avgr=0. ! <rho>
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) ! <c>
avgr=avgr/REAL(ny*nz) ! <rho>
avgc=avgc/REAL(ny*nz_all) ! <c>
avgr=avgr/REAL(ny*nz_all) ! <rho>
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) ! <c>
avgr=avgr/REAL(ny*nz) ! <rho>
avgc=avgc/REAL(ny*nz_all) ! <c>
avgr=avgr/REAL(ny*nz_all) ! <rho>
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

2
m_work.f90
View file

@ -13,7 +13,7 @@ module m_work
implicit none
! --- work arrays
real*4,allocatable :: tmp4(:,:,:)
real*8, allocatable :: tmp4(:,:,:)
real*8, allocatable :: wrk(:,:,:,:)

17
main.f90
View file

@ -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