Homec4science
Diffusion wireles-gpu (master) Browse

wireles-gpu/src/wireles_src.cuf
61fde96d0eeemaster

Last Change
Raw File

wireles_src.cuf
View Options

program wireles_src
use iso_c_binding
use precision
use dimen
use io
use cudafor
use poisson
use divergence
use cufft_m
use filter_gpu
use derivative
use wall_gpu
use forcing
use turbine_m
use convection
use sgs
use step_gpu
implicit none
!---------------------------------------------------------------------------
! iterator
integer*4 t
integer i,j,k
!---------------------------------------------------------------------------
! flow
real(fp),dimension(nx,ny):: zo_host,ustar_host,m_host
real(fp),dimension(nx,ny,nz2) :: u_host,v_host,w_host,&
rhsx_host,rhsy_host,rhsz_host
real(fp), dimension(nx,ny,nz2) :: lm_old_host,mm_old_host,qn_old_host,nn_old_host
real(fp),dimension(nx,ny,nz2) :: fx_host,fy_host,fz_host
real(fp),dimension(nx,ny),device :: zo,ustar,m
real(fp),dimension(nx,ny,nz2),device :: u,v,w,p,& ! Velocity & pressure
cx,cy,cz,& ! Convection term
txx,txy,txz,& ! Stress term
tyy,tyz,tzz,&
divtx,divty,divtz,& ! Divergence term
rhsx,rhsy,rhsz,& ! right hand side term
rhsx_f,rhsy_f,rhsz_f,&
dudx,dudy,dudz,& ! Velocity gradient
dvdx,dvdy,dvdz,&
dwdx,dwdy,dwdz,& ! Pressure gradient
dpdx,dpdy,dpdz,&
rhs_p,& ! Poisson eq rhs
lm_old,mm_old,&
qn_old,nn_old,cs2
!---------------------------------------------------------------------------
! turbine
real(fp), dimension(nx,ny,nz2),device :: fx,fy,fz
real(fp), dimension(3,181) :: alpha_cl_cd_host
real(fp), dimension(3,181),device :: alpha_cl_cd
!---------------------------------------------------------------------------
! io
real(fp) :: t01
real(fp) :: log_ustar, log_kemax, log_umax
real(fp), dimension(ta_nx,ta_ny,nzb) :: ta_u_host,ta_v_host,ta_w_host,&
ta_u2_host,ta_v2_host,ta_w2_host,&
ta_uw_host,ta_txz_host,ta_dudz_host
real(fp), dimension(ta_nx,ta_ny,nzb),device :: ta_u,ta_v,ta_w,&
ta_u2,ta_v2,ta_w2
!---------------------------------------------------------------------------
! fft plan
type(c_ptr) plan_batch(2),plan_2d(2)
integer istat
!---------------------------------------------------------------------------
! cuda variable
type(dim3) :: grid, tBlock, grid1,tBlock1,grid_2d,tBlock_2d
grid = dim3(nx/16,ny/16,nz2/2)
tBlock = dim3(16,16,2)
! Notice the difference of nzb vs nz2
! TODO: FINE TUNING THE BLOCK GRID CONFIGURATION
grid1 = dim3(nx/16,ny/16,nzb/2)
tBlock1 = dim3(16,16,2)
grid_2d = dim3(nx/16,ny/16,1)
tBlock_2d = dim3(16,16,1)
real(fp), device :: d_avg1,d_avg2
real(fp) :: d_avg_debug
integer :: device_id
!===========================================================================
! init
!===========================================================================
!---------------------------------------------------------------------------
! init mpi
!---------------------------------------------------------------------------
call mpi_init( ierr )
call mpi_comm_rank( mpi_comm_world, me, ierr )
call mpi_comm_size( mpi_comm_world, job_np, ierr )
nall=mpi_comm_world
if(me == 0) t01=mpi_wtime()
if (job_np /= nprocs) then
print *, 'error: job_np /= nprocs'
call exit()
endif
!---------------------------------------------------------------------------
! print start (important for debug)
!---------------------------------------------------------------------------
call print_start(me)
!---------------------------------------------------------------------------
! load variables
!---------------------------------------------------------------------------
call load_variable(u_host,'input/u.bin',me,nall)
call load_variable(v_host,'input/v.bin',me,nall)
call load_variable(w_host,'input/w.bin',me,nall)
! bc bot
if (me==0) then
w_host(:,:,2)=0._fp
end if
! bc top (part1)
if (me==nprocs-1) then
w_host(:,:,nzb+1)=0._fp
end if
! call update_uv(u,me,nall)
! call update_uv(v,me,nall)
! call update_w(w,me,nall)
! bc top (part2)
if (me==nprocs-1) then
u_host(:,:,nzb+1)=u_host(:,:,nzb)
v_host(:,:,nzb+1)=v_host(:,:,nzb)
endif
! Copy to device
u = u_host
v = v_host
w = w_host
if (resub_flag > 0 .or. sim_flag > 1) then
call load_variable(rhsx_host,'input/rhsx.bin',me,nall)
call load_variable(rhsy_host,'input/rhsy.bin',me,nall)
call load_variable(rhsz_host,'input/rhsz.bin',me,nall)
call load_variable(lm_old_host,'input/lm_old.bin',me,nall)
call load_variable(mm_old_host,'input/mm_old.bin',me,nall)
call load_variable(qn_old_host,'input/qn_old.bin',me,nall)
call load_variable(nn_old_host,'input/nn_old.bin',me,nall)
! Copy to device
rhsx = rhsx_host
rhsy = rhsy_host
rhsz = rhsz_host
lm_old = lm_old_host
mm_old = mm_old_host
qn_old = qn_old_host
nn_old = nn_old_host
endif
call load_surface(zo_host,'input/zo.bin')
zo = zo_host
!---------------------------------------------------------------------------
! load turbine aero data
!---------------------------------------------------------------------------
if(turb_flag > 0)then
open(unit=99,file='input/cl_cd_wire01.dat')
do k=1,181
read(99,*)alpha_cl_cd_host(1:3,k)
end do
close(99)
endif
! Copy to device
alpha_cl_cd = alpha_cl_cd_host
!---------------------------------------------------------------------------
! init io
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! init log
if(log_flag > 0) call open_log(me)
!---------------------------------------------------------------------------'
! init ta outputs
if (ta_flag > 0) then
if (ta_flag == 1) then
call initialize_ta1(ta_u,ta_v,ta_w,&
ta_u2,ta_v2,ta_w2)
call open_ta1(me)
! elseif (ta_flag == 2) then
! call initialize_ta2(ta_u,ta_v,ta_w,ta_u2,ta_v2,ta_w2,ta_uw,ta_txz,ta_dudz)
! call open_ta2(me)
end if
end if
print *, 'init ta outputs complete'
!---------------------------------------------------------------------------
! init ts outputs
if (ts_flag > 0) then
if (ts_flag == 1) then
call open_ts1(me)
! call open_ts2(me)
end if
end if
!---------------------------------------------------------------------------
! init inflow
if (sim_flag > 1) then
call open_inflow_m(me)
end if
!---------------------------------------------------------------------------
! init turb
if(turb_flag > 0) call open_turb(me)
!---------------------------------------------------------------------------
! init fft
! Watch out for the difference of row-first and column first indexing between fortran and cufft
if (fp==singlePrecision) then
call cufftPlanMany(plan_batch(1),2,(/ny,nx/),(/ny,nx/),1,nx*ny,(/ny,nx/2+1/),1,(nx/2+1)*ny,cufft_R2C,nzb)
call cufftPlanMany(plan_batch(2),2,(/ny,nx/),(/ny,nx/2+1/),1,(nx/2+1)*ny,(/ny,nx/),1,nx*ny,cufft_C2R,nzb)
call cufftPlan2d(plan_2d(1),ny,nx,CUFFT_R2C)
call cufftPlan2d(plan_2d(2),ny,nx,CUFFT_C2R)
else
call cufftPlanMany(plan_batch(1),2,(/ny,nx/),(/ny,nx/),1,nx*ny,(/ny,nx/2+1/),1,(nx/2+1)*ny,cufft_D2Z,nzb)
call cufftPlanMany(plan_batch(2),2,(/ny,nx/),(/ny,nx/2+1/),1,(nx/2+1)*ny,(/ny,nx/),1,nx*ny,cufft_Z2D,nzb)
call cufftPlan2d(plan_2d(1),ny,nx,CUFFT_D2Z)
call cufftPlan2d(plan_2d(2),ny,nx,CUFFT_Z2D)
end if
!---------------------------------------------------------------------------
! copy cpu data to gpu memory
device_id = me
!===========================================================================
! start time loop
!===========================================================================
do t=1,nsteps
!-----------------------------------------------------------------------
! save previous time's right-hand-sides for adams-bashforth
!-----------------------------------------------------------------------
rhsx_f = rhsx
rhsy_f = rhsy
rhsz_f = rhsz
!-----------------------------------------------------------------------
! derivatives
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! spectral deriv (and filter)
call ddxy_filter_gpu(u,dudx,dudy,t,plan_batch)
call ddxy_filter_gpu(v,dvdx,dvdy,t,plan_batch)
call ddxy_filter_gpu(w,dwdx,dwdy,t,plan_batch)
!-----------------------------------------------------------------------
! bc top uv
if (me==nprocs-1) then
!$cuf kernel do(2) <<<(nx/16,ny/16),(16,16)>>>
do j=1,ny
do i=1,nx
u(i,j,nzb+1)=u(i,j,nzb)
v(i,j,nzb+1)=v(i,j,nzb)
end do
end do
endif
!-----------------------------------------------------------------------
! spatial deriv
call ddz_gpu<<<grid, tBlock>>>(dudz,u,0) ! cell centre
call ddz_gpu<<<grid, tBlock>>>(dvdz,v,0) ! cell centre
call ddz_gpu<<<grid, tBlock>>>(dwdz,w,1) ! staggered
!-----------------------------------------------------------------------
! wall
!-----------------------------------------------------------------------
! we need to hide this wall modif to the other processor because
! dudz, dvdz change location (from w node to uv nodes) --> no update
if (me==0) then
! Porte-agel 2000 JFM correction
! TODO: RENAME THE VARIABLE: NOT A MEAN HERE
!$cuf kernel do(2) <<<(*,*),(16,16)>>>
do i = 1, nx
do j = 1, ny
d_avg1 = d_avg1 + dudz(i,j,3) ! Plane average
d_avg2 = d_avg2 + dvdz(i,j,3) ! Plane average
end do
end do
! ! Special treate ment for nzb == 1
! if (nzb > 1 .and. me == 0) then
! k_corr = 3
! elseif (nzb == 1 .and. me == 1) then
! k_corr = 2
! endif
call ddz_uv_corr_gpu<<<grid_2d, tBlock_2d>>>(dudz,d_avg1,3)
call ddz_uv_corr_gpu<<<grid_2d, tBlock_2d>>>(dvdz,d_avg2,3)
call wall_stress(txz,tyz,dudz,dvdz,ustar,m,u,v,zo,t,plan_2d)
end if
!-----------------------------------------------------------------------
! compute convec
!-----------------------------------------------------------------------
call convec_gpu<<<grid,tBlock>>>(cx,cy,cz,&
u,v,w,&
dudy,dudz,&
dvdx,dvdz,&
dwdx,dwdy )
!-----------------------------------------------------------------------
! compute sgs: LASD model
!-----------------------------------------------------------------------
call sgs_stag_gpu(u,v,w,&
dudx,dudy,dudz,&
dvdx,dvdy,dvdz,&
dwdx,dwdy,dwdz, &
LM_old,MM_old,QN_old,NN_old,cs2,&
txx,txy,txz,&
tyy,tyz,tzz,t,me,nall)
!-----------------------------------------------------------------------
! top bc
if (me==nprocs-1) then
!$cuf kernel do <<<*,*>>>
do j=1,ny
do i=1,nx
tzz (i,j,nzb+1)=tzz (i,j,nzb)
end do
end do
end if
!-----------------------------------------------------------------------
! Divergence
!-----------------------------------------------------------------------
call divstress_gpu(divtx,txx,txy,txz,1,t,plan_batch) ! w node
call divstress_gpu(divty,txy,tyy,tyz,1,t,plan_batch) ! w node
call divstress_gpu(divtz,txz,tyz,tzz,0,t,plan_batch) ! uv node
!-----------------------------------------------------------------------
! assemble matrix
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! rhs0
call rhs_assemble<<<grid, tBlock>>>(rhsx,cx,divtx )
call rhs_assemble<<<grid, tBlock>>>(rhsy,cy,divty )
call rhs_assemble<<<grid, tBlock>>>(rhsz,cz,divtz )
!-----------------------------------------------------------------------
! forcing
if (sim_flag < 3) call const_pres_grad<<<grid, tBlock>>>(rhsx )
if (turb_flag == 1) then
if (device_id == 0) then
call turbine_gpu(fx,fy,fz,u,v,w,alpha_cl_cd,t)
call turbine_force<<<grid, tBlock>>>(rhsx,fx )
call turbine_force<<<grid, tBlock>>>(rhsy,fy )
call turbine_force<<<grid, tBlock>>>(rhsz,fz )
end if
end if
!-----------------------------------------------------------------------
! Time stepping
!-----------------------------------------------------------------------
if (t == 1 .and. sim_flag < 2) then
call step<<<grid,tBlock>>>(u,rhsx,rhsx_f,1)
call step<<<grid,tBlock>>>(v,rhsy,rhsy_f,1)
call step<<<grid,tBlock>>>(w,rhsz,rhsz_f,1)
else
call step<<<grid,tBlock>>>(u,rhsx,rhsx_f,0)
call step<<<grid,tBlock>>>(v,rhsy,rhsy_f,0)
call step<<<grid,tBlock>>>(w,rhsz,rhsz_f,0)
end if
! bc bot w
if (me==0) then
!$cuf kernel do <<<*,*>>>
do j=1,ny
do i=1,nx
w (i,j,2)=0._fp
end do
end do
end if
! bc top w
if (me==nprocs-1) then
!$cuf kernel do <<<*,*>>>
do j=1,ny
do i=1,nx
w (i,j,nzb+1)=0._fp
end do
end do
end if
!-----------------------------------------------------------------------
! pressure correction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! get div
call ddxy(dudx,u,1,t,plan_batch)
call ddxy(dvdy,v,2,t,plan_batch)
call ddz_gpu<<<grid, tBlock>>>(dwdz,w,1)
call div_p<<<grid, tBlock>>>(rhs_p,dudx,dvdy,dwdz)
!-----------------------------------------------------------------------
! get pressure grandients (and press)
call press_stag_gpu(p,rhs_p,dpdx,dpdy,me,nall,t,plan_batch)
! ! ! call update_uv(p,me,nall)
call ddz_gpu<<<grid, tBlock>>>(dpdz,p,0)
! bc bot dpdz
if (me==0) then
!$cuf kernel do <<<*,*>>>
do j=1,ny
do i=1,nx
dpdz (i,j,2) = 0.0_fp
end do
end do
end if
! bc top dpdz
if (me==nprocs-1) then
!$cuf kernel do <<<*,*>>>
do j=1,ny
do i=1,nx
dpdz (i,j,nzb+1) = 0.0_fp
end do
end do
end if
!----------------------------------------------------------------------
! correct velocity
call velo_corr<<<grid, tBlock>>>(u,dpdx )
call velo_corr<<<grid, tBlock>>>(v,dpdy )
call velo_corr<<<grid, tBlock>>>(w,dpdz )
!-----------------------------------------------------------------------
! io
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! print log
if (mod(t,c_count*10)==0) then
u_host = u
v_host = v
w_host = w
fx_host = fx
fy_host = fy
fz_host = fz
ustar_host = ustar
call print_log(log_ustar,log_umax,u_host,v_host,w_host,ustar_host,t01,t,me,nall)
if(log_flag > 0) call write_log(log_ustar,log_umax,me)
if (log_umax>100) exit
if (ts_flag > 0 .and. t >= ts_tstart) then
print *,'Write instant field:'
call write_ts1(u_host,v_host,w_host,me)
end if
end if
! !-----------------------------------------------------------------------
! ! write inflow
if (sim_flag == 2) then
call write_inflow(u,v,w,t)
end if
! !-----------------------------------------------------------------------
! ! read (and force) inflow
if (sim_flag == 3) then
call read_inflow(u,v,w,t)
end if
! if (sim_flag > 2 .and. mod(t,2)== 0) then
! call read_inflow_m(u,v,w,me)
! end if
! !-----------------------------------------------------------------------
! ! compute time average outputs (just the sum at this stage)
if (ta_flag > 0 .and. t >= ta_tstart .and. mod(t,c_count) == 0) then
if (ta_flag==1) then
call compute_ta1<<<grid1,tBlock1>>>(u,v,w,&
ta_u,ta_v,ta_w,&
ta_u2,ta_v2,ta_w2 )
! elseif (ta_flag == 2) then
! call compute_ta2(u,v,w,txz,dudz,ta_u,ta_v,ta_w,ta_u2,ta_v2,ta_w2,ta_uw,ta_txz,ta_dudz,me)
end if
end if
! !-----------------------------------------------------------------------
! ! write and re-initalize ta_ (from the sum)
if (ta_flag > 0 .and. t >= ta_tstart .and. mod(t,p_count) == 0) then
print *,'Write average field:'
if (ta_flag==1) then
ta_u_host = ta_u
ta_v_host = ta_v
ta_w_host = ta_w
ta_u2_host = ta_u2
ta_v2_host = ta_v2
ta_w2_host = ta_w2
call write_ta1(ta_u_host,ta_v_host,ta_w_host,&
ta_u2_host,ta_v2_host,ta_w2_host,me)
call initialize_ta1(ta_u,ta_v,ta_w,&
ta_u2,ta_v2,ta_w2 )
! elseif (ta_flag == 2) then
! call write_ta2(ta_u,ta_v,ta_w,ta_u2,ta_v2,ta_w2,ta_uw,ta_txz,ta_dudz,me)
! call initialize_ta2(ta_u,ta_v,ta_w,ta_u2,ta_v2,ta_w2,ta_uw,ta_txz,ta_dudz)
end if
end if
! !-----------------------------------------------------------------------
! ! ts outputs
! !-----------------------------------------------------------------------
! ! turb outputs
! ! if (t >= ta_tstart .and. turb_flag > 0 .and. me == 0) then
! ! ta_thrust = ta_thrust + thrust
! ! ta_power = ta_power + power
! ! if (mod(t,p_count) == 0) then
! ! print *, 'ta_thrust: ', ta_thrust/p_count
! ! print *, 'ta_power: ', ta_power/p_count
! ! ta_thrust = 0._fp
! ! ta_power = 0._fp
! ! endif
! ! endif
enddo
! ===============================================================================
! end time loop
!===============================================================================
!---------------------------------------------------------------------------
! close log
!---------------------------------------------------------------------------
if(log_flag > 0) call close_log(me)
! !---------------------------------------------------------------------------
! ! close inflow
! !---------------------------------------------------------------------------
if (sim_flag > 1) then
call close_inflow_m(me)
end if
! !---------------------------------------------------------------------------
! ! close and recompose ta outputs
! !---------------------------------------------------------------------------
if (ta_flag > 0) then
if (ta_flag==1) then
call close_ta1(me)
call recompose_ta1(me)
! elseif (ta_flag == 2) then
! call close_ta2(me)
! call recompose_ta2(me)
end if
end if
! !---------------------------------------------------------------------------
! ! close and recompose ts outputs
! !---------------------------------------------------------------------------
if (ts_flag > 0) then
if (ts_flag==1) then
call close_ts1(me)
! call close_ts2(me)
call recompose_ts1(me)
end if
end if
!---------------------------------------------------------------------------
! close turb
!---------------------------------------------------------------------------
! if(turb_flag > 0) call close_turb(me)
!---------------------------------------------------------------------------
! save variables
!---------------------------------------------------------------------------
u_host = u
v_host = v
w_host = w
lm_old_host = lm_old
mm_old_host = mm_old
qn_old_host = qn_old
nn_old_host = nn_old
rhsx_host = rhsx
rhsy_host = rhsy
rhsz_host = rhsz
call save_variable(u_host,'input/u.bin',me,nall)
call save_variable(v_host,'input/v.bin',me,nall)
call save_variable(w_host,'input/w.bin',me,nall)
call save_variable(lm_old_host,'input/lm_old.bin',me,nall)
call save_variable(mm_old_host,'input/mm_old.bin',me,nall)
call save_variable(qn_old_host,'input/qn_old.bin',me,nall)
call save_variable(nn_old_host,'input/nn_old.bin',me,nall)
call save_variable(rhsx_host,'input/rhsx.bin',me,nall)
call save_variable(rhsy_host,'input/rhsy.bin',me,nall)
call save_variable(rhsz_host,'input/rhsz.bin',me,nall)
!---------------------------------------------------------------------------
! save warmup data for precursor and main
!---------------------------------------------------------------------------
if (me == 0 .and. sim_flag < 2) then
if (resub_flag == 0 ) call system('mkdir ../init_data')
call system('cp input/* ../init_data')
end if
!---------------------------------------------------------------------------
! print end
!---------------------------------------------------------------------------
call print_end(t01,me)
!---------------------------------------------------------------------------
! mpi finalize
!---------------------------------------------------------------------------
999 call mpi_finalize(ierr)
end program wireles_src
c4science · Help