Ignore:
Timestamp:
Mar 8, 2013 11:54:10 PM (11 years ago)
Author:
raasch
Message:

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:


resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:


bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/fft_xy.f90

    r1107 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! further openACC statements added, CUDA branch completely runs on GPU
     23! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
     24! used total domain)
    2325!
    2426! Former revisions:
     
    213215          total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
    214216          total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
    215           CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (ny+1)*nz )
    216           CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (ny+1)*nz )
    217           CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nx+1)*nz )
    218           CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nx+1)*nz )
     217          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
     218          CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
     219          CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
     220          CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
    219221#else
    220222          message_string = 'no system-specific fft-call available'
     
    261263
    262264       CHARACTER (LEN=*) ::  direction
    263        INTEGER ::  i, ishape(1), j, k, m
     265       INTEGER ::  i, ishape(1), j, k
    264266
    265267       LOGICAL ::  forward_fft
     
    273275       REAL, DIMENSION(6*(nx+1)) ::  work2
    274276#elif defined( __cuda_fft )
    275        REAL(dpk), DEVICE, DIMENSION(:), ALLOCATABLE    ::  cuda_a_device
    276        COMPLEX(dpk), DEVICE, DIMENSION(:), ALLOCATABLE ::  cuda_b_device
    277        COMPLEX(dpk), DIMENSION(:), ALLOCATABLE         ::  cuda_host
     277       !$acc declare create( ar_tmp )
     278       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
    278279#endif
    279280       REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ar
     
    502503#elif defined( __cuda_fft )
    503504
    504           ALLOCATE( cuda_a_device(0:total_points_x_transpo-1) )
    505           ALLOCATE( cuda_b_device(0:((nx+1)/2+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) - 1) )
    506           ALLOCATE( cuda_host(0:((nx+1)/2+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) - 1) )
    507 
    508           m = 0
    509 
    510505          IF ( forward_fft )  THEN
    511506
    512              cuda_a_device = ar(0:total_points_x_transpo-1,nys_x,nzb_x)
    513 
    514              CALL CUFFTEXECD2Z( plan_xf, cuda_a_device, cuda_b_device )
    515              cuda_host = cuda_b_device
    516 
     507             !$acc data present( ar )
     508             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
     509
     510             !$acc kernels
     511             !$acc loop
    517512             DO  k = nzb_x, nzt_x
    518513                DO  j = nys_x, nyn_x
    519514
     515                   !$acc loop vector( 32 )
    520516                   DO  i = 0, (nx+1)/2
    521                       ar(i,j,k)      = REAL( cuda_host(m+i) )  * dnx
    522                    ENDDO
    523 
     517                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
     518                   ENDDO
     519
     520                   !$acc loop vector( 32 )
    524521                   DO  i = 1, (nx+1)/2 - 1
    525                       ar(nx+1-i,j,k) = AIMAG( cuda_host(m+i) ) * dnx
    526                    ENDDO
    527 
    528                    m = m + (nx+1)/2 + 1
    529 
    530                 ENDDO
    531              ENDDO
    532 
    533           ELSE
    534 
     522                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
     523                   ENDDO
     524
     525                ENDDO
     526             ENDDO
     527             !$acc end kernels
     528             !$acc end data
     529
     530          ELSE
     531
     532             !$acc data present( ar )
     533             !$acc kernels
     534             !$acc loop
    535535             DO  k = nzb_x, nzt_x
    536536                DO  j = nys_x, nyn_x
    537537
    538                    cuda_host(m) = CMPLX( ar(0,j,k), 0.0 )
    539 
     538                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
     539
     540                   !$acc loop vector( 32 )
    540541                   DO  i = 1, (nx+1)/2 - 1
    541                       cuda_host(m+i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
    542                    ENDDO
    543                    cuda_host(m+(nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
    544 
    545                    m = m + (nx+1)/2 + 1
    546 
    547                 ENDDO
    548              ENDDO
    549 
    550              cuda_b_device = cuda_host
    551              CALL CUFFTEXECZ2D( plan_xi, cuda_b_device, cuda_a_device )
    552 
    553              ar(0:total_points_x_transpo-1,nys_x,nzb_x) = cuda_a_device
    554 
    555           ENDIF
    556 
    557           DEALLOCATE( cuda_a_device, cuda_b_device, cuda_host )
     542                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
     543                   ENDDO
     544                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
     545
     546                ENDDO
     547             ENDDO
     548             !$acc end kernels
     549
     550             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
     551             !$acc end data
     552
     553          ENDIF
    558554
    559555#else
     
    775771
    776772       CHARACTER (LEN=*) ::  direction
    777        INTEGER ::  i, j, jshape(1), k, m
     773       INTEGER ::  i, j, jshape(1), k
    778774
    779775       LOGICAL ::  forward_fft
     
    787783       REAL, DIMENSION(6*(ny+1)) ::  work2
    788784#elif defined( __cuda_fft )
    789        REAL(dpk), DEVICE, DIMENSION(:), ALLOCATABLE    ::  cuda_a_device
    790        COMPLEX(dpk), DEVICE, DIMENSION(:), ALLOCATABLE ::  cuda_b_device
    791        COMPLEX(dpk), DIMENSION(:), ALLOCATABLE         ::  cuda_host
     785       !$acc declare create( ar_tmp )
     786       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
    792787#endif
    793788       REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar
     
    10131008#elif defined( __cuda_fft )
    10141009
    1015           ALLOCATE( cuda_a_device(0:total_points_y_transpo-1) )
    1016           ALLOCATE( cuda_b_device(0:((ny+1)/2+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) - 1) )
    1017           ALLOCATE( cuda_host(0:((ny+1)/2+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) - 1) )
    1018 
    1019           m = 0
    1020 
    10211010          IF ( forward_fft )  THEN
    10221011
    1023              cuda_a_device = ar(0:total_points_y_transpo-1,nxl_y,nzb_y)
    1024 
    1025              CALL CUFFTEXECD2Z( plan_yf, cuda_a_device, cuda_b_device )
    1026              cuda_host = cuda_b_device
    1027 
     1012             !$acc data present( ar )
     1013             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
     1014
     1015             !$acc kernels
     1016             !$acc loop
    10281017             DO  k = nzb_y, nzt_y
    10291018                DO  i = nxl_y, nxr_y
    10301019
     1020                   !$acc loop vector( 32 )
    10311021                   DO  j = 0, (ny+1)/2
    1032                       ar(j,i,k)      = REAL( cuda_host(m+j) )  * dny
    1033                    ENDDO
    1034 
     1022                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
     1023                   ENDDO
     1024
     1025                   !$acc loop vector( 32 )
    10351026                   DO  j = 1, (ny+1)/2 - 1
    1036                       ar(ny+1-j,i,k) = AIMAG( cuda_host(m+j) ) * dny
    1037                    ENDDO
    1038 
    1039                    m = m + (ny+1)/2 + 1
    1040 
    1041                 ENDDO
    1042              ENDDO
    1043 
    1044           ELSE
    1045 
     1027                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
     1028                   ENDDO
     1029
     1030                ENDDO
     1031             ENDDO
     1032             !$acc end kernels
     1033             !$acc end data
     1034
     1035          ELSE
     1036
     1037             !$acc data present( ar )
     1038             !$acc kernels
     1039             !$acc loop
    10461040             DO  k = nzb_y, nzt_y
    10471041                DO  i = nxl_y, nxr_y
    10481042
    1049                    cuda_host(m) = CMPLX( ar(0,i,k), 0.0 )
    1050 
     1043                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
     1044
     1045                   !$acc loop vector( 32 )
    10511046                   DO  j = 1, (ny+1)/2 - 1
    1052                       cuda_host(m+j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
    1053                    ENDDO
    1054                    cuda_host(m+(ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
    1055 
    1056                    m = m + (ny+1)/2 + 1
    1057 
    1058                 ENDDO
    1059              ENDDO
    1060 
    1061              cuda_b_device = cuda_host
    1062              CALL CUFFTEXECZ2D( plan_yi, cuda_b_device, cuda_a_device )
    1063 
    1064              ar(0:total_points_y_transpo-1,nxl_y,nzb_y) = cuda_a_device
    1065 
    1066           ENDIF
    1067 
    1068           DEALLOCATE( cuda_a_device, cuda_b_device, cuda_host )
     1047                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
     1048                   ENDDO
     1049                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
     1050
     1051                ENDDO
     1052             ENDDO
     1053             !$acc end kernels
     1054
     1055             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
     1056             !$acc end data
     1057
     1058          ENDIF
    10691059
    10701060#else
Note: See TracChangeset for help on using the changeset viewer.