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/dissip
Files:
1 edited
1 moved

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/dissip/dissip_gcm.F90

    r933 r953  
    11MODULE dissip_gcm_mod 
    22  USE icosa 
    3  
     3  USE abort_mod 
    44  PRIVATE 
    55 
    … …  
    4848  USE icosa 
    4949  IMPLICIT NONE   
    50     CALL allocate_field(f_due_diss1,field_u,type_real,llm) 
    51     CALL allocate_field(f_due_diss2,field_u,type_real,llm) 
     50    CALL allocate_field(f_due_diss1,field_u,type_real,llm,ondevice=.TRUE.) 
     51    CALL allocate_field(f_due_diss2,field_u,type_real,llm,ondevice=.TRUE.) 
    5252    CALL allocate_field(f_dtheta_diss,field_t,type_real,llm) 
    53     CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm) 
     53    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm,ondevice=.TRUE.) 
    5454    ALLOCATE(tau_graddiv(llm)) 
    5555    ALLOCATE(tau_gradrot(llm))     
    5656    ALLOCATE(tau_divgrad(llm)) 
     57    !$acc enter data create(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async 
    5758  END SUBROUTINE allocate_dissip 
    5859  
    … …  
    6667  USE transfert_omp_mod 
    6768  USE omp_para 
     69  USE abort_mod 
    6870  IMPLICIT NONE 
    6971   
    … …  
    98100      IF (is_master) PRINT *, 'No Rayleigh friction' 
    99101   CASE('dcmip2_schaer_noshear') 
     102      CALL abort_acc("rayleigh_friction_type /= 'none'") 
    100103      rayleigh_friction_type=1 
    101104      rayleigh_shear=0 
    102105      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1' 
    103106   CASE('dcmip2_schaer_shear') 
     107      CALL abort_acc("rayleigh_friction_type /= 'none'") 
    104108      rayleigh_shear=1 
    105109      rayleigh_friction_type=2 
    106110      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' 
    107111   CASE('giant_liu_schneider')  
     112      CALL abort_acc("rayleigh_friction_type /= 'none'") 
    108113      rayleigh_friction_type=99 
    109114      IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010' 
    … …  
    130135 
    131136    CALL allocate_dissip 
    132     CALL allocate_field(f_u,field_u,type_real) 
    133     CALL allocate_field(f_du,field_u,type_real) 
    134     CALL allocate_field(f_theta,field_t,type_real) 
    135     CALL allocate_field(f_dtheta,field_t,type_real) 
     137    CALL allocate_field(f_u,field_u,type_real,ondevice=.TRUE.) 
     138    CALL allocate_field(f_du,field_u,type_real,ondevice=.TRUE.) 
     139    CALL allocate_field(f_theta,field_t,type_real,ondevice=.TRUE.) 
     140    CALL allocate_field(f_dtheta,field_t,type_real,ondevice=.TRUE.) 
    136141     
    137142    CALL init_message(f_due_diss1,req_e1_vect,req_due) 
    … …  
    173178      CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
    174179 
     180      ! This cannot be ported on GPU due to compiler limitations 
    175181      DO j=jj_begin,jj_end 
    176182        DO i=ii_begin,ii_end 
    … …  
    192198      dumax=0 
    193199      DO iter=1,nitergdiv 
     200        CALL update_device_field(f_u) 
    194201        CALL transfert_request(f_u,req_e1_vect) 
     202          
    195203        DO ind=1,ndomain 
    196204          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    … …  
    200208          du=f_du(ind) 
    201209          CALL compute_gradiv_inplace(u,1,1) 
     210          ! This should be ported on GPU but we are running into compiler issues... 
     211          !$acc update host(u(:)) wait 
    202212          du=u 
    203213        ENDDO 
    204214      ENDDO 
    205215 
     216      CALL update_device_field(f_du) 
    206217      CALL transfert_request(f_du,req_e1_vect) 
    207        
     218      CALL update_host_field(f_du) 
     219 
    208220      DO ind=1,ndomain 
    209221        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    210222        CALL swap_dimensions(ind) 
    211223        CALL swap_geometry(ind) 
    212         u=f_u(ind) 
    213224        du=f_du(ind) 
    214225          
     226        ! Not ported on GPU because all the other kernels cannot be ported 
    215227        DO j=jj_begin,jj_end 
    216228          DO i=ii_begin,ii_end 
    … …  
    240252        u=f_u(ind) 
    241253        du=f_du(ind) 
     254        ! This should be ported on GPU but we are running into compiler issues... 
    242255        u=du/dumax 
    243256      ENDDO 
    … …  
    265278      CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
    266279  
     280      ! This cannot be ported on GPU due to compiler limitations 
    267281       DO j=jj_begin,jj_end 
    268282        DO i=ii_begin,ii_end 
    … …  
    280294!$OMP BARRIER 
    281295 
    282  
    283296    DO it=1,20 
    284297  
    285298      dumax=0 
    286299      DO iter=1,nitergrot 
     300        CALL update_device_field(f_u) 
    287301        CALL transfert_request(f_u,req_e1_vect) 
    288302        DO ind=1,ndomain 
    … …  
    293307          du=f_du(ind) 
    294308          CALL compute_gradrot_inplace(u,1,1) 
     309          ! This should be ported on GPU but we are running into compiler issues... 
     310          !$acc update host(u(:)) wait 
    295311          du=u 
    296312        ENDDO 
    297313      ENDDO 
    298  
     314       
     315      CALL update_device_field(f_du) 
    299316      CALL transfert_request(f_du,req_e1_vect) 
     317      CALL update_host_field(f_du) 
    300318       
    301319      DO ind=1,ndomain 
    … …  
    303321        CALL swap_dimensions(ind) 
    304322        CALL swap_geometry(ind) 
    305         u=f_u(ind) 
    306323        du=f_du(ind) 
    307324         
     325        ! Not ported on GPU because all the other kernels cannot be ported 
    308326        DO j=jj_begin,jj_end 
    309327          DO i=ii_begin,ii_end 
    … …  
    334352        u=f_u(ind) 
    335353        du=f_du(ind) 
     354        ! This should be ported on GPU but we are running into compiler issues... 
    336355        u=du/dumax 
    337356      ENDDO 
    … …  
    358377      CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 
    359378 
     379      ! This cannot be ported on GPU due to compiler limitations 
    360380      DO j=jj_begin,jj_end 
    361381        DO i=ii_begin,ii_end 
    … …  
    373393      dthetamax=0 
    374394      DO iter=1,niterdivgrad 
     395        CALL update_device_field(f_theta) 
    375396        CALL transfert_request(f_theta,req_i1) 
    376397        DO ind=1,ndomain 
    … …  
    381402          dtheta=f_dtheta(ind) 
    382403          CALL compute_divgrad_inplace(theta,1,1) 
     404          ! This should be ported on GPU but we are running into compiler issues... 
     405          !$acc update host(theta(:)) wait 
    383406          dtheta=theta 
    384407        ENDDO 
    385408      ENDDO 
    386409 
     410      CALL update_device_field(f_dtheta) 
    387411      CALL transfert_request(f_dtheta,req_i1) 
     412      CALL update_host_field(f_dtheta) 
    388413       
    389414      DO ind=1,ndomain 
    … …  
    391416        CALL swap_dimensions(ind) 
    392417        CALL swap_geometry(ind) 
    393         theta=f_theta(ind) 
    394418        dtheta=f_dtheta(ind) 
    395419         
     420        ! Not ported on GPU because all the other kernels cannot be ported 
    396421        DO j=jj_begin,jj_end 
    397422          DO i=ii_begin,ii_end 
    … …  
    421446        theta=f_theta(ind) 
    422447        dtheta=f_dtheta(ind) 
     448        ! This should be ported on GPU but we are running into compiler issues... 
    423449        theta=dtheta/dthetamax 
    424450      ENDDO 
    … …  
    489515    ENDIF 
    490516 
     517    !$acc update device(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async 
     518 
    491519  END SUBROUTINE init_dissip  
    492520  
    … …  
    501529  USE time_mod 
    502530  USE omp_para 
     531  USE abort_mod 
    503532  IMPLICIT NONE 
    504533    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:) 
    … …  
    534563      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind) 
    535564 
     565      !$acc parallel loop collapse(2) present(due(:,:), dtheta_rhodz(:,:,:), due_diss1(:,:), due_diss2(:,:), dtheta_rhodz_diss(:,:), tau_graddiv(:), tau_gradrot(:), tau_divgrad(:)) async 
    536566      DO l=ll_begin,ll_end 
    537567!DIR$ SIMD 
    … …  
    545575        ENDDO 
    546576      ENDDO 
     577      !$acc end parallel loop 
    547578 
    548579!      dtheta_rhodz=0 
    … …  
    550581 
    551582      IF(rayleigh_friction_type>0) THEN 
     583        CALL abort_acc("dissip/rayleigh_friction_type>0") 
    552584       IF(rayleigh_friction_type<99) THEN 
    553585         phi=f_geopot(ind) 
    … …  
    562594       ELSE 
    563595         ue=f_ue(ind) 
     596            !$acc parallel loop present(ue(:,:), due(:,:), lat_e(:)) 
    564597            DO ij=ij_begin,ij_end 
    565598              nn = ij+u_right 
    … …  
    577610              ENDIF 
    578611            ENDDO 
     612            !$acc end parallel loop 
    579613       ENDIF 
    580614      END IF 
    581615   END DO 
    582  
    583616   CALL trace_end("dissip") 
    584617 
    585 !   CALL write_dissip_tendencies 
     618   !CALL write_dissip_tendencies 
    586619!$OMP BARRIER 
    587620 
    … …  
    653686      ue=f_ue(ind) 
    654687      due=f_due(ind)  
     688      !$acc parallel loop present(ue(:,:), due(:,:)) async 
    655689      DO  l = ll_begin, ll_end 
     690        !$acc loop 
    656691!DIR$ SIMD 
    657692        DO ij=ij_begin,ij_end 
    … …  
    664699 
    665700    DO it=1,nitergdiv 
    666          
    667701      CALL send_message(f_due,req_due) 
    668702      CALL wait_message(req_due) 
    669          
     703 
    670704      DO ind=1,ndomain 
    671705        IF (.NOT. assigned_domain(ind)) CYCLE 
    672706        CALL swap_dimensions(ind) 
    673707        CALL swap_geometry(ind) 
    674         due=f_due(ind)  
     708        due=f_due(ind) 
    675709        CALL compute_gradiv_inplace(due,ll_begin,ll_end) 
    676710      ENDDO 
    … …  
    702736      ue=f_ue(ind) 
    703737      due=f_due(ind)  
     738      !$acc parallel loop present(ue(:,:), due(:,:)) async 
    704739      DO  l = ll_begin, ll_end 
     740         !$acc loop 
    705741!DIR$ SIMD 
    706742        DO ij=ij_begin,ij_end 
    … …  
    713749 
    714750    DO it=1,nitergrot 
    715          
    716751      CALL send_message(f_due,req_due) 
    717752      CALL wait_message(req_due) 
    … …  
    740775    REAL(rstd),POINTER    :: theta(:,:) 
    741776    REAL(rstd),POINTER    :: dtheta(:,:) 
    742     INTEGER :: ind 
     777    INTEGER :: ind, l, ij 
    743778    INTEGER :: it 
    744779 
    … …  
    751786      theta=f_theta(ind) 
    752787      dtheta=f_dtheta(ind)  
    753       dtheta=theta 
     788      ! Replace Fortran 90 construct dtheta=theta because it confuses PGI acc kernels... 
     789      !$acc parallel loop collapse(2) present(theta(:,:), dtheta(:,:)) async 
     790      DO l=ll_begin,ll_end 
     791!DIR$ SIMD 
     792        DO ij=1,iim*jjm 
     793          dtheta(ij,l)=theta(ij,l) 
     794        ENDDO 
     795      ENDDO 
     796      !$acc end parallel loop 
    754797    ENDDO 
    755798 
    756799    DO it=1,niterdivgrad 
    757          
    758800      CALL transfert_request(f_dtheta,req_i1) 
    759801         
    … …  
    797839      theta_rhodz=f_theta_rhodz(ind) 
    798840      dtheta_rhodz=f_dtheta_rhodz(ind)  
     841      !$acc parallel loop present(mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:)) async 
    799842      DO  l = ll_begin, ll_end 
     843        !$acc loop 
    800844!DIR$ SIMD 
    801845        DO ij=ij_begin,ij_end 
    … …  
    806850 
    807851    DO it=1,niterdivgrad 
    808          
    809852      CALL send_message(f_dtheta_rhodz,req_dtheta) 
    810853      CALL wait_message(req_dtheta) 
    … …  
    826869      mass=f_mass(ind) 
    827870     
     871      !$acc parallel loop collapse(2) present(mass(:,:), dtheta_rhodz(:,:)) async 
    828872      DO  l = ll_begin, ll_end 
    829873!DIR$ SIMD 
    … …  
    832876        ENDDO 
    833877      ENDDO 
     878      !$acc end parallel loop 
    834879    ENDDO 
    835880 
    … …  
    837882    CALL trace_end("divgrad") 
    838883 
    839   END SUBROUTINE divgrad_theta_rhodz   
    840    
     884  END SUBROUTINE divgrad_theta_rhodz 
     885 
    841886  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle) 
    842   USE icosa 
    843   IMPLICIT NONE 
    844887    INTEGER,INTENT(IN)     :: llb 
    845888    INTEGER,INTENT(IN)     :: lle 
     889    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm) 
    846890    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm) 
    847     REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm) 
    848      
     891 
    849892    gradivu_e = ue 
    850     CALL compute_gradiv_inplace(gradivu_e, llb, lle) 
     893    CALL compute_gradiv_inplace(gradivu_e,llb,lle) 
     894 
    851895  END SUBROUTINE compute_gradiv 
    852896   
    853897  SUBROUTINE compute_gradiv_inplace(ue_gradivu_e,llb,lle) 
    854   USE icosa 
    855   IMPLICIT NONE 
     898    USE geometry, ONLY : Ai, ne, le, de 
    856899    INTEGER,INTENT(IN)     :: llb 
    857900    INTEGER,INTENT(IN)     :: lle 
    … …  
    861904    INTEGER :: ij,l 
    862905 
     906    ! ue and gradivu_e second dimension is not always llm so use the bounds explicitly 
     907    !$acc data present( ue_gradivu_e(:,llb:lle), Ai(:), ne(:,:), le(:), de(:)) create(divu_i(:,:)) async 
     908 
     909    !$acc parallel loop collapse(2) async 
    863910    DO l=llb,lle 
    864911      !DIR$ SIMD 
    … …  
    872919      ENDDO 
    873920    ENDDO 
     921    !$acc end parallel loop 
    874922     
     923    !$acc parallel loop collapse(2) async 
    875924    DO l=llb,lle 
    876925      !DIR$ SIMD 
    … …  
    881930      ENDDO 
    882931    ENDDO 
    883  
     932    !$acc end parallel loop 
     933 
     934    !$acc parallel loop collapse(2) async 
    884935    DO l=llb,lle 
    885936      !DIR$ SIMD 
    … …  
    889940          ue_gradivu_e(ij+u_ldown,l)=-ue_gradivu_e(ij+u_ldown,l)*cgraddiv 
    890941      ENDDO 
    891     ENDDO    
     942    ENDDO 
     943    !$acc end parallel loop 
     944    !$acc end data    
    892945  END SUBROUTINE compute_gradiv_inplace 
    893    
     946 
    894947  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle) 
    895   USE icosa 
    896   IMPLICIT NONE 
    897948    INTEGER,INTENT(IN)     :: llb 
    898949    INTEGER,INTENT(IN)     :: lle 
    899950    REAL(rstd),INTENT(IN) :: theta(iim*jjm,1:lle) 
    900951    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,1:lle) 
    901      
     952 
    902953    divgrad_i = theta 
    903     CALL compute_divgrad_inplace(divgrad_i, llb, lle) 
    904   END SUBROUTINE 
     954    CALL compute_divgrad_inplace(divgrad_i,llb,lle) 
     955  END SUBROUTINE compute_divgrad 
    905956   
    906957  SUBROUTINE compute_divgrad_inplace(theta_divgrad_i,llb,lle) 
    907   USE icosa 
    908   IMPLICIT NONE 
     958    USE geometry, ONLY : Ai, ne, le, de 
    909959    INTEGER,INTENT(IN)     :: llb 
    910960    INTEGER,INTENT(IN)     :: lle 
    … …  
    912962    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle) 
    913963 
    914     INTEGER :: ij,l     
     964    INTEGER :: ij,l 
     965 
     966    ! theta and divgrad_i second dimension is not always llm so use the bounds explicitly 
     967    !$acc data present(theta_divgrad_i(:,llb:lle), Ai(:), de(:), ne(:,:), le(:)) create(grad_e(:,:)) async 
     968        
     969    !$acc parallel loop collapse(2) async 
    915970    DO l=llb,lle 
    916971      !DIR$ SIMD 
    … …  
    920975          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta_divgrad_i(ij,l)+ne(ij+t_ldown,rup)*theta_divgrad_i(ij+t_ldown,l) ) 
    921976      ENDDO 
    922     ENDDO     
     977    ENDDO 
     978    !$acc end parallel loop 
    923979     
     980     
     981    !$acc parallel loop collapse(2) async 
    924982    DO l=llb,lle 
    925983      !DIR$ SIMD 
    … …  
    933991      ENDDO 
    934992    ENDDO 
     993    !$acc end parallel loop 
    935994    
     995    !$acc parallel loop collapse(2) async 
    936996    DO l=llb,lle 
    937997      DO ij=ij_begin,ij_end 
    … …  
    939999      ENDDO 
    9401000    ENDDO 
     1001    !$acc end parallel loop 
     1002    !$acc end data 
    9411003  END SUBROUTINE compute_divgrad_inplace 
    9421004 
    … …  
    9461008    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,lle) 
    9471009    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,lle) 
    948      
     1010 
    9491011    gradrot_e = ue 
    9501012    CALL compute_gradrot_inplace(gradrot_e,llb,lle) 
    951   END SUBROUTINE 
    952      
     1013  END SUBROUTINE compute_gradrot 
     1014 
    9531015  SUBROUTINE compute_gradrot_inplace(ue_gradrot_e,llb,lle) 
    954   USE icosa 
    955   IMPLICIT NONE 
     1016    USE geometry, ONLY : Av, ne, le, de 
    9561017    INTEGER,INTENT(IN)     :: llb 
    9571018    INTEGER,INTENT(IN)     :: lle 
    … …  
    9611022    INTEGER :: ij,l 
    9621023 
     1024    ! ue and gradrot_e second dimension is not always llm so use the bounds explicitly 
     1025    ! gradrot_e should be copyout but using copy instead allows to compare the output 
     1026    ! more easily as the code sometimes uses unintialed values 
     1027    !$acc data present(ue_gradrot_e(:,llb:lle), Av(:), ne(:,:), de(:), le(:)) create(rot_v(:,:)) async 
     1028      
     1029    !$acc parallel loop collapse(2) async 
    9631030    DO l=llb,lle 
    9641031!DIR$ SIMD 
    … …  
    9721039      ENDDO 
    9731040    ENDDO                               
    974    
     1041    !$acc end parallel loop 
     1042   
     1043    !$acc parallel loop collapse(2) async 
    9751044    DO l=llb,lle 
    9761045!DIR$ SIMD 
    … …  
    9811050      ENDDO 
    9821051    ENDDO 
    983  
     1052    !$acc end parallel loop 
     1053 
     1054    !$acc parallel loop collapse(2) async 
    9841055    DO l=llb,lle 
    9851056!DIR$ SIMD 
    … …  
    9891060          ue_gradrot_e(ij+u_ldown,l)=-ue_gradrot_e(ij+u_ldown,l)*cgradrot         
    9901061      ENDDO 
    991     ENDDO   
     1062    ENDDO 
     1063    !$acc end parallel loop 
     1064    !$acc end data    
    9921065  END SUBROUTINE compute_gradrot_inplace 
    9931066 
    9941067 
    9951068END MODULE dissip_gcm_mod 
    996         

  • codes/icosagcm/trunk/src/dissip/sponge.f90

    r548 r953  
    2424  USE omp_para 
    2525  USE mpipara, ONLY: is_mpi_master 
     26  USE abort_mod 
    2627  IMPLICIT NONE 
    2728    INTEGER :: l 
    … …  
    4344      RETURN 
    4445    ENDIF 
     46 
     47    IF (iflag_sponge > 0) THEN 
     48       CALL abort_acc("iflag_sponge > 0") 
     49    END IF 
    4550 
    4651!$OMP MASTER 
Note: See TracChangeset for help on using the changeset viewer.