diff --git a/init_velocity.f90 b/init_velocity.f90 index ff2f690..f9e4f5f 100644 --- a/init_velocity.f90 +++ b/init_velocity.f90 @@ -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 diff --git a/m_force.f90 b/m_force.f90 index da0f8bf..f14e993 100644 --- a/m_force.f90 +++ b/m_force.f90 @@ -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