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/time
Files:
3 moved

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/time/euler_scheme.F90

    r933 r953  
    11MODULE euler_scheme_mod 
    22  USE field_mod 
     3  USE abort_mod 
    34  IMPLICIT NONE 
    45  PRIVATE 
    … …  
    4546 
    4647       IF(with_dps) THEN ! update ps/mass 
     48          CALL abort_acc("Euler_scheme/with_dps") 
    4749          IF(caldyn_eta==eta_mass) THEN ! update ps 
    4850             ps=f_ps(ind) ; dps=f_dps(ind) ;               
    … …  
    7173       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 
    7274 
    73        DO l=ll_begin,ll_end 
     75       CALL compute_euler_scheme(u, du, theta_rhodz, dtheta_rhodz) 
     76    ENDDO 
     77 
     78    CALL trace_end("Euler_scheme") 
     79 
     80    CONTAINS 
     81 
     82    SUBROUTINE compute_euler_scheme(u, du, theta_rhodz, dtheta_rhodz) 
     83       REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm) 
     84       REAL(rstd),INTENT(IN)    :: du(iim*3*jjm,llm) 
     85       REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn) 
     86       REAL(rstd),INTENT(IN)    :: dtheta_rhodz(iim*jjm,llm,nqdyn) 
     87 
     88       !$acc data present(theta_rhodz(:,:,:), u(:,:),du(:,:), dtheta_rhodz(:,:,:)) async 
     89 
     90       !$acc parallel loop async 
     91       DO l=ll_begin,ll_end 
     92          !$acc loop 
    7493          !DIR$ SIMD 
    7594          DO ij=ij_begin,ij_end 
    … …  
    8099          ENDDO 
    81100       ENDDO 
    82     ENDDO 
    83  
    84     CALL trace_end("Euler_scheme")   
     101        
     102       !$acc end data 
     103    END SUBROUTINE compute_euler_scheme 
    85104 
    86105  END SUBROUTINE Euler_scheme 
    … …  
    97116    INTEGER :: l,ij 
    98117 
     118    !$acc data present(hflux(:,:), wflux(:,:), hfluxt(:,:), wfluxt(:,:)) async 
     119 
    99120    IF(fluxt_zero) THEN 
    100121 
    101122       fluxt_zero=.FALSE. 
    102  
    103        DO l=ll_begin,ll_end 
     123       !$acc parallel loop async 
     124       DO l=ll_begin,ll_end 
     125          !$acc loop 
    104126          !DIR$ SIMD 
    105127          DO ij=ij_begin_ext,ij_end_ext 
    … …  
    111133 
    112134       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
     135          !$acc parallel loop async 
    113136          DO l=ll_begin,ll_endp1 
     137             !$acc loop 
    114138             !DIR$ SIMD 
    115139             DO ij=ij_begin,ij_end 
    … …  
    120144 
    121145    ELSE 
    122  
    123        DO l=ll_begin,ll_end 
     146       !$acc parallel loop async 
     147       DO l=ll_begin,ll_end 
     148          !$acc loop 
    124149          !DIR$ SIMD 
    125150          DO ij=ij_begin_ext,ij_end_ext 
    … …  
    131156 
    132157       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
     158          !$acc parallel loop async 
    133159          DO l=ll_begin,ll_endp1 
     160             !$acc loop 
    134161             !DIR$ SIMD 
    135162             DO ij=ij_begin,ij_end 
    … …  
    138165          ENDDO 
    139166       END IF 
    140  
    141167    END IF 
    142  
     168    !$acc end data 
    143169  END SUBROUTINE accumulate_fluxes 
    144170 
    145171  SUBROUTINE legacy_to_DEC(f_ps, f_u) 
     172    USE icosa 
     173    USE disvert_mod 
     174    USE omp_para 
     175    USE trace 
     176    TYPE(t_field),POINTER :: f_ps(:), f_u(:) 
     177    REAL(rstd), POINTER :: ps(:), u(:,:) 
     178    INTEGER :: ind,ij,l 
     179 
     180    CALL trace_start("legacy_to_DEC") 
     181 
     182    DO ind=1,ndomain 
     183       IF (.NOT. assigned_domain(ind)) CYCLE 
     184       CALL swap_dimensions(ind) 
     185       CALL swap_geometry(ind) 
     186 
     187       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps 
     188          ps=f_ps(ind) 
     189          !$acc parallel loop async default(present) 
     190          !DIR$ SIMD 
     191          DO ij=ij_begin,ij_end 
     192             ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass 
     193          ENDDO 
     194       END IF 
     195 
     196       u=f_u(ind) 
     197       !$acc parallel loop async default(present) 
     198       DO l=ll_begin,ll_end 
     199          !$acc loop 
     200          !DIR$ SIMD 
     201          DO ij=ij_begin,ij_end 
     202             u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right) 
     203             u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup) 
     204             u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown) 
     205          ENDDO 
     206       ENDDO 
     207    ENDDO 
     208 
     209    CALL trace_end("legacy_to_DEC")   
     210 
     211  END SUBROUTINE Legacy_to_DEC 
     212 
     213  SUBROUTINE DEC_to_legacy(f_ps, f_u) 
    146214    USE icosa 
    147215    USE disvert_mod 
    … …  
    158226       CALL swap_geometry(ind) 
    159227 
    160        IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps 
    161           ps=f_ps(ind) 
    162           !DIR$ SIMD 
    163           DO ij=ij_begin,ij_end 
    164              ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass 
    165           ENDDO 
    166        END IF 
    167  
    168        u=f_u(ind) 
    169        DO l=ll_begin,ll_end 
    170           !DIR$ SIMD 
    171           DO ij=ij_begin,ij_end 
    172              u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right) 
    173              u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup) 
    174              u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown) 
    175           ENDDO 
    176        ENDDO 
    177     ENDDO 
    178  
    179     CALL trace_end("legacy_to_DEC")   
    180   END SUBROUTINE Legacy_to_DEC 
    181  
    182   SUBROUTINE DEC_to_legacy(f_ps, f_u) 
    183     USE icosa 
    184     USE disvert_mod 
    185     USE omp_para 
    186     USE trace 
    187     TYPE(t_field),POINTER :: f_ps(:), f_u(:) 
    188     REAL(rstd), POINTER :: ps(:), u(:,:) 
    189     INTEGER :: ind,ij,l 
    190     CALL trace_start("legacy_to_DEC") 
    191  
    192     DO ind=1,ndomain 
    193        IF (.NOT. assigned_domain(ind)) CYCLE 
    194        CALL swap_dimensions(ind) 
    195        CALL swap_geometry(ind) 
    196  
    197228       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN 
    198229          ps=f_ps(ind) 
     230          !$acc parallel loop async default(present) 
    199231          !DIR$ SIMD 
    200232          DO ij=ij_begin,ij_end 
    … …  
    204236        
    205237       u=f_u(ind) 
    206        DO l=ll_begin,ll_end 
     238       !$acc parallel loop async default(present) 
     239       DO l=ll_begin,ll_end 
     240          !$acc loop 
    207241          !DIR$ SIMD 
    208242          DO ij=ij_begin,ij_end 

  • codes/icosagcm/trunk/src/time/hevi_scheme.F90

    r933 r953  
    44  USE field_mod 
    55  USE euler_scheme_mod 
     6  USE abort_mod 
    67  IMPLICIT NONE 
    78  PRIVATE 
    … …  
    1314 
    1415CONTAINS 
    15  
    1616  SUBROUTINE set_coefs_ark23(dt) 
    1717    ! ARK2 scheme by Giraldo, Kelly, Constantinescu 2013 
    … …  
    6060 
    6161    CALL legacy_to_DEC(f_ps, f_u) 
     62 
    6263    DO j=1,nb_stage 
    6364       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), & 
    … …  
    7677          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
    7778          CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, wj(j), fluxt_zero(ind)) 
     79 
    7880       END DO 
    7981       ! update model state 
    … …  
    8890          CALL update_3D(cjl(l,j), f_u, f_du_fast(:,l)) 
    8991          IF(.NOT. hydrostatic) THEN 
     92             CALL abort_acc("HEVI_scheme/!hydrostatic") 
    9093             CALL update_3D(bjl(l,j), f_W, f_dW_slow(:,l)) 
    9194             CALL update_3D(cjl(l,j), f_W, f_dW_fast(:,l)) 
    … …  
    9396             CALL update_3D(cjl(l,j), f_geopot, f_dPhi_fast(:,l)) 
    9497          END IF 
     98 
    9599!$OMP BARRIER        
    96100       END DO 
    … …  
    142146    INTENT(IN) :: dy 
    143147    INTEGER :: l 
     148    !$acc kernels async default(present) 
    144149    DO l=ll_begin,ll_end 
    145150       y(:,l)=y(:,l)+w*dy(:,l) 
    146151    ENDDO 
     152    !$acc end kernels 
    147153  END SUBROUTINE compute_update_3D 
    148154   
    … …  
    164170    INTENT(INOUT) :: y 
    165171    INTENT(IN) :: dy 
     172    !$acc kernels async default(present) 
    166173    y(:)=y(:)+w*dy(:) 
     174    !$acc end kernels 
    167175  END SUBROUTINE compute_update_2D 
    168176   

  • codes/icosagcm/trunk/src/time/timeloop_gcm.F90

    r933 r953  
    9494    END SELECT 
    9595     
     96    IF (scheme_family /= hevi) THEN 
     97       CALL abort_acc("scheme_family /= hevi") 
     98    END IF 
     99 
    96100    ! Time-independant orography 
    97101    CALL allocate_field(f_phis,field_t,type_real,name='phis') 
    … …  
    103107    CALL allocate_field(f_u,field_u,type_real,llm,name='u') 
    104108    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') 
    105     CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') 
     109    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') ! used only if .not. hydrostatic 
    106110    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q') 
    107111    ! Mass fluxes 
    108     CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes 
    109     CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time 
     112    CALL allocate_field(f_hflux,field_u,type_real,llm, ondevice=.TRUE.)    ! instantaneous mass fluxes 
     113    CALL allocate_field(f_hfluxt,field_u,type_real,llm,ondevice=.TRUE.)   ! mass "fluxes" accumulated in time 
    110114    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes 
    111     CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 
     115    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt',ondevice=.TRUE.) 
    112116     
    113117    SELECT CASE(scheme_family) 
    … …  
    125129    CASE(hevi) 
    126130       ! Trends 
    127        CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow') 
    128        CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow') 
    129        CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast') 
    130        CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow') 
    131        CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast') 
     131       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow', ondevice=.TRUE.) 
     132       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow', ondevice=.TRUE.) 
     133       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast', ondevice=.TRUE.) 
     134       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow', ondevice=.TRUE.) 
     135       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast', ondevice=.TRUE.) 
    132136       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow') 
    133137       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast') 
    … …  
    172176   
    173177  SUBROUTINE timeloop 
     178    USE abort_mod 
    174179    USE dissip_gcm_mod 
    175180    USE sponge_mod 
    … …  
    211216       rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind) 
    212217       IF(caldyn_eta==eta_mass) THEN 
    213           CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 
     218          CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.FALSE.) ! save rhodz for transport scheme before dynamics update ps 
    214219       ELSE 
    215220          DO l=ll_begin,ll_end 
    … …  
    244249    CALL SYSTEM_CLOCK(start_clock, rate_clock) 
    245250    !$OMP END MASTER    
     251    call update_device_field(f_ps) 
     252    call update_device_field(f_mass) 
     253    CALL update_device_field(f_theta_rhodz) 
     254    CALL update_device_field(f_u) 
     255    CALL update_device_field(f_q) 
     256    CALL update_device_field(f_geopot) 
     257    CALL update_device_field(f_wflux) 
     258    CALL update_device_field(f_rhodz) 
     259 
    246260 
    247261    DO it=itau0+1,itau0+itaumax 
    … …  
    263277          CALL wait_message(req_mass0) 
    264278          CALL send_message(f_theta_rhodz,req_theta_rhodz0)  
    265           CALL wait_message(req_theta_rhodz0)  
     279          CALL wait_message(req_theta_rhodz0) 
    266280          CALL send_message(f_u,req_u0) 
    267281          CALL wait_message(req_u0) 
    … …  
    281295       SELECT CASE(scheme_family) 
    282296       CASE(explicit) 
     297          CALL abort_acc("explicit_scheme") 
    283298          CALL explicit_scheme(it, fluxt_zero) 
    284299       CASE(hevi) 
    … …  
    298313                CALL swap_geometry(ind) 
    299314                mass=f_mass(ind); ps=f_ps(ind); 
    300                 CALL compute_rhodz(.TRUE., ps, mass) 
     315                CALL compute_rhodz(.TRUE., ps, mass, ondevice=.TRUE.) 
    301316             END DO 
    302317          ENDIF 
    … …  
    311326          CALL euler_scheme(.FALSE.)  ! update only u, theta 
    312327          IF (iflag_sponge > 0) THEN 
     328             CALL abort_acc("iflag_sponge>0") 
    313329             CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz) 
    314330             CALL euler_scheme(.FALSE.)  ! update only u, theta 
    … …  
    321337       END IF 
    322338       CALL exit_profile(id_dissip) 
    323         
     339 
    324340       CALL enter_profile(id_adv) 
    325341       IF(MOD(it,itau_adv)==0) THEN 
    … …  
    329345          ! At this point advect_tracer has obtained the halos of u and rhodz, 
    330346          ! needed for correct computation of kinetic energy 
     347          IF(diagflux_on) CALL abort_acc("diagflux_on") 
    331348          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_buf_i, f_hfluxt) 
    332349 
    … …  
    340357             END DO 
    341358          ENDIF 
     359          IF(positive_theta) CALL abort_acc("positive_theta") 
    342360          IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q) 
    343361       END IF 
    344362       CALL exit_profile(id_adv) 
    345         
     363 
    346364       CALL enter_profile(id_diags) 
    347365!       IF (MOD(it,itau_physics)==0) THEN 
    … …  
    360378 
    361379       IF (MOD(it,itau_check_conserv)==0) THEN 
     380          CALL update_host_field(f_ps) 
     381          CALL update_host_field(f_theta_rhodz) 
     382          CALL update_host_field(f_u) 
     383          CALL update_host_field(f_dps) 
     384          CALL update_host_field(f_q) 
    362385          CALL check_conserve_detailed(it, AAM_dyn, & 
    363386               f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
    … …  
    367390       IF (mod(it,itau_out)==0 ) THEN 
    368391          CALL transfert_request(f_u,req_e1_vect) 
     392          CALL update_host_field(f_ps)               
     393          CALL update_host_field(f_mass) 
     394          CALL update_host_field(f_theta_rhodz) 
     395          CALL update_host_field(f_geopot) 
     396          CALL update_host_field(f_u) 
     397          CALL update_host_field(f_q) 
    369398          CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q) 
    370399       ENDIF 
    … …  
    374403    END DO 
    375404     
     405    CALL update_host_field(f_ps) 
     406    CALL update_host_field(f_theta_rhodz) 
     407    CALL update_host_field(f_u) 
     408    CALL update_host_field(f_q) 
     409    CALL update_host_field(f_geopot) 
     410 
    376411!    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  
    377412    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)  
    378      
     413 
     414    CALL update_host_field(f_dps)     
    379415    CALL check_conserve_detailed(it, AAM_dyn, & 
    380416         f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
Note: See TracChangeset for help on using the changeset viewer.