Ignore:
Timestamp:
07/15/19 12:29:31 (5 years ago)
Author:
adurocher
Message:

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )

Location:
codes/icosagcm/trunk/src/transport
Files:
1 edited
1 moved

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/transport/advect.F90

    r899 r953  
    4949  !======================================================================================= 
    5050 
    51   SUBROUTINE compute_gradq3d(qi,sqrt_leng,gradq3d,xyz_i,xyz_v) 
     51  SUBROUTINE compute_gradq3d(qi,sqrt_leng,gradq3d,xyz_i, xyz_v) 
    5252    USE trace 
    5353    USE omp_para 
    … …  
    6868    REAL(rstd) :: x1,x2,x3 
    6969    REAL(rstd) :: dq(3) 
    70  
     70     
    7171    CALL trace_start("compute_gradq3d1") 
    7272 
    … …  
    8686!    END DO 
    8787 
    88      DO l = ll_begin,ll_end  
     88     !$acc data create(gradtri(:,:,:), arr(:), ar(:)) present(sqrt_leng(:), xyz_i(:,:), xyz_v(:,:), qi(:,:), gradq3d(:,:,:)) async 
     89 
     90     !$acc parallel loop collapse(2) async private(A, dq) 
     91     DO l = ll_begin,ll_end 
    8992!DIR$ SIMD 
    9093      DO ij=ij_begin_ext,ij_end_ext 
    … …  
    151154         
    152155      ENDDO 
     156     ENDDO 
    153157       
     158     !$acc parallel loop collapse(2) async private(A, dq) 
     159     DO l = ll_begin,ll_end 
    154160      DO ij=ij_begin_ext,ij_end_ext 
    155161 
    … …  
    219225 
    220226!DIR$ SIMD 
     227    !$acc parallel loop async 
    221228    DO ij=ij_begin,ij_end 
    222229       ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50 
    … …  
    225232    CALL trace_start2("compute_gradq3d2") 
    226233       
     234    !$acc parallel loop collapse(3) async 
    227235    DO k=1,3 
    228236      DO l =ll_begin,ll_end 
    … …  
    239247 
    240248    !============================================================================================= LIMITING  
     249    !$acc parallel loop collapse(2) async 
    241250    DO l =ll_begin,ll_end 
    242251!DIR$ SIMD 
    … …  
    251260             minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), & 
    252261                  qi(ij+t_rdown,l),qi(ij+t_ldown,l)) 
    253              alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) ) 
    254              alphamx = max(alphamx,0.0) 
    255              alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l)) 
    256              alphami = max(alphami,0.0)  
     262             IF ((maxq_c - qi(ij,l)) /= 0.0) THEN 
     263               alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) ) 
     264               alphamx = max(alphamx,0.0) 
     265             ELSE 
     266               alphamx = 0.0 
     267             ENDIF 
     268             IF ((minq_c - qi(ij,l)) /= 0.0) THEN 
     269               alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l)) 
     270               alphami = max(alphami,0.0) 
     271             ELSE 
     272               alphami = 0.0 
     273             ENDIF 
    257274             alpha   = min(alphamx,alphami,1.0) 
    258275!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 
    … …  
    264281 
    265282  CALL trace_end("compute_gradq3d3") 
     283 
     284  !$acc end data 
    266285   
    267286  CONTAINS 
    … …  
    329348 
    330349  ! Backward trajectories, for use with Miura approach 
    331   SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc) 
     350  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau,cc) 
     351    USE geometry, ONLY : xyz_e, de, wee, le 
    332352    USE trace 
    333353    USE omp_para 
    … …  
    345365 
    346366    ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement 
    347      
     367 
     368    !$acc data present(ue(:,:), cc(:,:,:), normal(:,:), tangent(:,:), xyz_e(:,:), de(:), wee(:,:,:), le(:)) async 
     369 
    348370    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau 
     371    !$acc parallel loop private(up_e, v_e) collapse(2) gang vector async 
    349372    DO l = ll_begin,ll_end 
    350373!DIR$ SIMD 
    … …  
    397420       ENDDO 
    398421    END DO 
    399  
     422    !$acc end data 
    400423    CALL trace_end("compute_backward_traj") 
    401424 
    … …  
    404427  ! Horizontal transport (S. Dubey, T. Dubos) 
    405428  ! Slope-limited van Leer approach with hexagons 
    406   SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass,qi,qfluxt) 
     429  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass, qi, qfluxt) 
    407430    USE trace 
    408431    USE omp_para 
     432    USE abort_mod 
     433    USE geometry, only : Ai, xyz_i 
    409434    IMPLICIT NONE 
    410435    LOGICAL, INTENT(IN)       :: update_mass, diagflux_on 
    … …  
    415440    REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm) 
    416441    REAL(rstd), INTENT(INOUT) :: qfluxt(3*iim*jjm,MERGE(llm,1,diagflux_on)) ! time-integrated tracer flux 
    417  
    418     REAL(rstd) :: dq,dmass,qe, newmass 
     442! metrics terms 
     443     
     444    REAL(rstd) :: dq,dmass,qe,newmass 
    419445    REAL(rstd) :: qflux(3*iim*jjm,llm) 
    420     INTEGER :: ij,l 
     446    INTEGER :: ij,l,ij_tmp 
     447 
     448    IF(diagflux_on) CALL abort_acc("compute_advect_horiz : diagflux_on") 
    421449 
    422450    CALL trace_start("compute_advect_horiz") 
    423451#include "../kernels/advect_horiz.k90" 
    424452    CALL trace_end("compute_advect_horiz") 
     453 
    425454  END SUBROUTINE compute_advect_horiz 
    426455 

  • codes/icosagcm/trunk/src/transport/advect_tracer.F90

    r933 r953  
    3232    INTEGER :: ind 
    3333 
    34     CALL allocate_field(f_normal,field_u,type_real,3, name='normal') 
    35     CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent') 
    36     CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d') 
    37     CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc') 
    38     CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng') 
    39     CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw') 
    40     CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw') 
    41     CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq') 
    42     CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq') 
     34    CALL allocate_field(f_normal,field_u,type_real,3, name='normal',ondevice=.TRUE.) 
     35    CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent',ondevice=.TRUE.) 
     36    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d',ondevice=.TRUE.) 
     37    CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc',ondevice=.TRUE.) 
     38    CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng',ondevice=.TRUE.) 
     39    CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw',ondevice=.TRUE.) 
     40    CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw',ondevice=.TRUE.) 
     41    CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq',ondevice=.TRUE.) 
     42    CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq',ondevice=.TRUE.) 
    4343     
    4444    DO ind=1,ndomain 
    … …  
    5252    END DO 
    5353 
     54    CALL update_device_field(f_tangent) 
     55    CALL update_device_field(f_normal) 
     56    CALL update_device_field(f_sqrt_leng) 
     57 
    5458  END SUBROUTINE init_advect_tracer 
    5559 
    … …  
    6064    USE write_field_mod 
    6165    USE tracer_mod 
     66    USE abort_mod 
    6267    TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux 
    6368    TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux 
    … …  
    138143       END DO 
    139144 
    140        CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc)  
    141  
     145       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
    142146    END DO 
    143147 
    … …  
    174178 
    175179          IF(frac>0.) THEN ! accumulate mass, mass flux and tracer mass 
     180             CALL abort_acc("frac>0.") 
    176181             qmasst  = f_qmasst(ind)  
     182             !$acc kernels default(present) async 
    177183             qmasst(:,ll_begin:ll_end,k) = qmasst(:,ll_begin:ll_end,k) + & 
    178184                  frac*rhodz(:,ll_begin:ll_end)*q(:,ll_begin:ll_end,k) 
     185             !$acc end kernels 
    179186             IF(k==nq_last) THEN 
    180187                masst  = f_masst(ind) 
    181188                massfluxt  = f_massfluxt(ind) 
     189                !$acc kernels default(present) async 
    182190                masst(:,ll_begin:ll_end) = masst(:,ll_begin:ll_end)+frac*rhodz(:,ll_begin:ll_end) 
    183191                massfluxt(:,ll_begin:ll_end) = massfluxt(:,ll_begin:ll_end)+hfluxt(:,ll_begin:ll_end) 
     192                !$acc end kernels 
    184193             END IF 
    185194          END IF 
    … …  
    189198      ENDIF 
    190199    END DO  
    191      
     200 
    192201    ! 1/2 vertical transport 
    193202!!$OMP BARRIER 
    … …  
    208217         IF (advection_scheme(k)==advect_vanleer) CALL vlz(k==nq_last, 0.5,wfluxt,rhodz, q(:,:,k),0, dzqw, adzqw, dzq, wq) 
    209218       END DO 
    210  
    211219    END DO 
    212220 
    … …  
    251259    ! finite difference of q 
    252260 
     261     !$acc data present(q(:,:), mass(:,:), wfluxt(:,:), dzqw(:,:), adzqw(:,:), dzq(:,:), wq(:,:)) async 
     262 
     263     !$acc parallel loop async collapse(2) 
    253264     DO l=ll_beginp1,ll_end 
    254265!DIR$ SIMD 
    … …  
    265276    ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l 
    266277 
     278     !$acc parallel loop async collapse(2) 
    267279     DO l=ll_beginp1,ll_endm1 
    268280!DIR$ SIMD 
    … …  
    281293    ! 0 slope in top and bottom layers 
    282294    IF (is_omp_first_level) THEN 
     295      !$acc parallel loop async 
    283296      DO ij=ijb,ije 
    284297           dzq(ij,1)=0. 
    … …  
    287300       
    288301    IF (is_omp_last_level) THEN 
     302      !$acc parallel loop async 
    289303      DO ij=ijb,ije 
    290304          dzq(ij,llm)=0. 
    … …  
    297311    ! sigw = fraction of mass that leaves level l/l+1 
    298312    ! then amount of q leaving level l/l+1 = wq = w * qq 
     313     !$acc parallel loop collapse(2) async 
    299314     DO l=ll_beginp1,ll_end 
    300315!DIR$ SIMD 
    … …  
    313328    ! wq = 0 at top and bottom 
    314329    IF (is_omp_first_level) THEN 
     330      !$acc parallel loop async 
    315331       DO ij=ijb,ije 
    316332            wq(ij,1)=0. 
    … …  
    319335     
    320336    IF (is_omp_last_level) THEN 
     337      !$acc parallel loop async 
    321338      DO ij=ijb,ije 
    322339            wq(ij,llm+1)=0. 
    … …  
    329346 
    330347    ! update q, mass is updated only after all q's have been updated 
     348    !$acc parallel loop collapse(2) async 
    331349    DO l=ll_begin,ll_end 
    332350!DIR$ SIMD 
    … …  
    337355       ENDDO 
    338356    END DO 
    339  
     357    !$acc end data 
    340358    CALL trace_end("vlz") 
    341359 
Note: See TracChangeset for help on using the changeset viewer.