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

FC-DFS: Fully Controlled Deterministic Forcing Scheme

Browse source
This commit is contained in:
ignis 2014-04-28 14:29:51 +09:00
parent 0a380a543f
commit abdc5f5fdc
2 changed files with 116 additions and 7 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
init_velocity.f90
View file

@ -208,7 +208,6 @@ subroutine init_velocity
! Ref. Flow, Turbulence and Combustion 73: 133-167, 2004
! by L. Guichard, J. Reveillon and R. Hauguel
! ------------------------------------------------------------------------
allocate(model_e_spec(kmax)) ; model_e_spec=zip
! open(unit=404,file="model_spectrum.dat")
! write(404,*) 'VARIABLES = "X","Energy(k)"' !
@ -218,6 +217,7 @@ subroutine init_velocity
! write(502,*) 'I=',kmax,'F=POINT'
if (isp_type.eq.4) then !ksj
model_e_spec=zip
k0=(SQRT(two*PI))/t_length
model_A=16.0*(SQRT(two/PI))*(t_u_prime**2.0)/k0
do i=1,kmax
@ -238,7 +238,6 @@ subroutine init_velocity
! e_spec_log01,e_spec_log02,e_spec_log03
! write(502,'(i5,e15.5)') i,model_e_spec(i)
enddo
endif
! close(404)
@ -290,6 +289,7 @@ subroutine init_velocity
! 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
!---------------------------------------------------------J. Kwon. 201403

119
m_force.f90
View file

@ -1,7 +1,8 @@
module m_force
use m_parameters, only : kfmax, famp, force_type
use m_parameters
implicit none
integer*4 :: n_forced_nodes, n_forced_nodes_total
@ -18,6 +19,8 @@ contains
implicit none
integer :: i, j, k, n, n_shell
real*8 :: dissi_spec(kmax)
real*8 :: eps_v_temp,eta_temp,uvar_temp,sctmp_temp,x_length_temp,u_prime_temp,tau_eta_temp
! if flow is not forced, return
if (flow_type .ne. 1) return
@ -26,9 +29,14 @@ contains
select case (force_type)
case (1)
case (1:2)
if (force_type == 1) then
! Machiels forcing (see article in PRL #79(18) p.3411)
write(out,*) "Forcing #1: Machiels forcing - setting up"
write(out,*) "Forcing #1: Machiels forcing - setting up"
else
! Fully Controlled Deterministic Forcing Scheme (FC-DFS)
write(out,*) "Forcing #2: FC-DFS - setting up"
end if
call flush(out)
! find out how many nodes are we forcing and book them
@ -52,8 +60,8 @@ contains
MPI_INTEGER4,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
! writing out the # of forced nodes
write(out,*) 'Number of forced nodes for this process:',n_forced_nodes
if (myid.eq.0) write(out,*) ' total number:',n_forced_nodes_total
write(out,'(a40,3x,i10)') 'Number of forced nodes for this process:',n_forced_nodes
if (myid.eq.0) write(out,'(a40,3x,i10)') ' total number:',n_forced_nodes_total
call flush(out)
! allocating arrays for the coordinates of the forced nodes
@ -77,6 +85,17 @@ contains
end do
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
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))
enddo
endif
!!$ ! writing out the nodes
!!$ do n = 1,n_forced_nodes
!!$ write(out,"(3i4)") ifn(n),jfn(n),kfn(n)
@ -116,6 +135,12 @@ contains
integer :: i, j, k, n_shell, n
real*8 :: fac, fac2
real*8 :: kk
real*8 :: mean_eps,tau_f
real*8 :: kfmax8
real*8 :: F_k(kmax),ff_k(kmax),eps_f(kmax)
real*8 :: alpha,sum_F_Em,sum_F_Es
real*8 :: energy_model
select case (force_type)
@ -172,6 +197,90 @@ contains
end do
!------------------------------------------------------------------------
! case (2)
! New forcing scheme added by ksj and edited by J. Kwon, 2014
! Ref. Flow, Turbulence and Combustion 73: 133-167, 2004
! by L. Guichard, J. Reveillon and R. Hauguel
! ------------------------------------------------------------------------
case (2)
! Fully Controlled Deterministic Forcing Scheme (FC-DFS)
e_spec = zip
e_spec1 = zip
hits = 0
hits1 = 0
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
! assembling the total energy in each shell
do k = 1,nz
do j = 1,ny
do i = 1,nx
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .gt. 0 .and. n_shell .le. kmax) then
fac2 = fac * (fields(i,j,k,1)**2 + fields(i,j,k,2)**2 + fields(i,j,k,3)**2)
if (akx(i).eq.0.d0) fac2 = fac2 * 0.5d0
e_spec1(n_shell) = e_spec1(n_shell) + fac2
end if
end do
end do
end do
! reducing the number of hits and energy to two arrays on master node
count = kmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
count = kmax
call MPI_BCAST(e_spec,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
! getting the total energy by integrating the spectrum
energy = sum(e_spec(1:kmax))
mean_eps = zip
tau_f = zip
eps_f = zip
F_k = zip
ff_k = zip
energy_model = sum(model_e_spec(1:kmax))
do i = 1,kmax
eps_f(i) = e_spec(i) * real(i**2,8) * two * nu
end do
mean_eps = sum(eps_f(1:kmax))
tau_f=c_f*sqrt(nu/mean_eps)
kfmax8=real(kfmax,8)
do i=1,kfmax
kk=real(i,8)
F_k(i)=(1.-exp(-((kk-kfmax8)**2.)/(delta_k**2.)))/(1.-exp(-(kfmax8**2.)/(delta_k**2.)))
enddo
sum_F_Em=zip
sum_F_Es=zip
do i=1,kfmax
sum_F_Em=sum_F_Em+F_k(i)*model_e_spec(i)
sum_F_Es=sum_F_Es+F_k(i)*e_spec(i)
enddo
alpha=(energy_model-energy+tau_f*mean_eps+sum_F_Es)/sum_F_Em
do i=1,kfmax
ff_k(i)=F_k(i)/(two*e_spec(i))*(alpha*model_e_spec(i)-e_spec(i))
enddo
do n = 1,n_forced_nodes
n_shell = k_shell(n)
i = ifn(n)
j = jfn(n)
k = kfn(n)
wrk(i,j,k,1) = wrk(i,j,k,1) + ff_k(n_shell)/tau_f * fields(i,j,k,1)
wrk(i,j,k,2) = wrk(i,j,k,2) + ff_k(n_shell)/tau_f * fields(i,j,k,2)
wrk(i,j,k,3) = wrk(i,j,k,3) + ff_k(n_shell)/tau_f * fields(i,j,k,3)
end do
case default
write(out,*) "WRONG FORCE_TYPE:",force_type