Ignore:
Timestamp:
Nov 8, 2013 3:18:40 PM (8 years ago)
Author:
raasch
Message:

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:


openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:


bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

File:
1 edited

Legend:

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

    r1093 r1257  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! openacc porting
     23! bugfix for calculation of advective timestep in case of vertically stretched
     24! grids
    2325!
    2426! Former revisions:
     
    9395    IMPLICIT NONE
    9496
    95     INTEGER ::  i, j, k, u_max_cfl_ijk(3), v_max_cfl_ijk(3)
    96 
    97     REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy, dt_plant_canopy_l,  &
    98              dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,      &
    99              dt_u, dt_v, dt_w, u_max_cfl, value, v_max_cfl
     97    INTEGER ::  i, j, k
     98
     99    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy, dt_plant_canopy_l,     &
     100             dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,         &
     101             dt_u, dt_u_l, dt_v, dt_v_l, dt_w, dt_w_l, u_gtrans_l, u_max_l,   &
     102             u_min_l, value, v_gtrans_l, v_max_l, v_min_l, w_max_l, w_min_l
    100103
    101104    REAL, DIMENSION(2)         ::  uv_gtrans, uv_gtrans_l
     105    REAL, DIMENSION(3)         ::  reduce, reduce_l
    102106    REAL, DIMENSION(nzb+1:nzt) ::  dxyz2_min
    103107
     
    127131!
    128132!--       Averaging over the entire model domain.
    129           uv_gtrans_l = 0.0
     133          u_gtrans_l = 0.0
     134          v_gtrans_l = 0.0
     135          !$acc parallel present( u, v )
    130136          DO  i = nxl, nxr
    131137             DO  j = nys, nyn
    132138                DO  k = nzb+1, nzt
    133                    uv_gtrans_l(1) = uv_gtrans_l(1) + u(k,j,i)
    134                    uv_gtrans_l(2) = uv_gtrans_l(2) + v(k,j,i)
     139                   u_gtrans_l = u_gtrans_l + u(k,j,i)
     140                   v_gtrans_l = v_gtrans_l + v(k,j,i)
    135141                ENDDO
    136142             ENDDO
    137143          ENDDO
    138           uv_gtrans_l = uv_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
     144          !$acc end parallel
     145          uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
     146          uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
    139147#if defined( __parallel )
    140148          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     
    151159
    152160!
    153 !-- Determine the maxima of the velocity components.
     161!-- Determine the maxima of the velocity components, including their
     162!-- grid index positions.
     163#if defined( __openacc )
     164    IF ( dt_fixed )  THEN  ! otherwise do it further below for better cache usage
     165       u_max_l = -999999.9
     166       u_min_l =  999999.9
     167       v_max_l = -999999.9
     168       v_min_l =  999999.9
     169       w_max_l = -999999.9
     170       w_min_l =  999999.9
     171       !$acc parallel present( u, v, w )
     172       DO  i = nxl, nxr
     173          DO  j = nys, nyn
     174             DO  k = nzb+1, nzt
     175                u_max_l = MAX( u_max_l, u(k,j,i) )
     176                u_min_l = MIN( u_min_l, u(k,j,i) )
     177                v_max_l = MAX( v_max_l, v(k,j,i) )
     178                v_min_l = MIN( v_min_l, v(k,j,i) )
     179                w_max_l = MAX( w_max_l, w(k,j,i) )
     180                w_min_l = MIN( w_min_l, w(k,j,i) )
     181             ENDDO
     182          ENDDO
     183       ENDDO
     184       !$acc end parallel
     185#if defined( __parallel )
     186       reduce_l(1) = u_max_l
     187       reduce_l(2) = v_max_l
     188       reduce_l(3) = w_max_l
     189       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     190       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
     191       u_max = reduce(1)
     192       v_max = reduce(2)
     193       w_max = reduce(3)
     194       reduce_l(1) = u_min_l
     195       reduce_l(2) = v_min_l
     196       reduce_l(3) = w_min_l
     197       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     198       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
     199       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
     200       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
     201       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
     202#else
     203       IF ( ABS( u_min_l ) > u_max_l )  THEN
     204          u_max = u_min_l
     205       ELSE
     206          u_max = u_max_l
     207       ENDIF
     208       IF ( ABS( v_min_l ) > v_max_l )  THEN
     209          v_max = v_min_l
     210       ELSE
     211          v_max = v_max_l
     212       ENDIF
     213       IF ( ABS( w_min_l ) > w_max_l )  THEN
     214          w_max = w_min_l
     215       ELSE
     216          w_max = w_max_l
     217       ENDIF
     218#endif
     219    ENDIF
     220#else
    154221    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0, &
    155222                         u_max, u_max_ijk )
     
    158225    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0, &
    159226                         w_max, w_max_ijk )
    160 
    161 !
    162 !-- In case of Galilei transformation, the horizontal velocity maxima have
    163 !-- to be calculated from the transformed horizontal velocities
    164     IF ( galilei_transformation )  THEN
    165        CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'absoff', &
    166                             u_gtrans, u_max_cfl, u_max_cfl_ijk )
    167        CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'absoff', &
    168                             v_gtrans, v_max_cfl, v_max_cfl_ijk )
    169     ELSE
    170        u_max_cfl = u_max
    171        v_max_cfl = v_max
    172        u_max_cfl_ijk = u_max_ijk
    173        v_max_cfl_ijk = v_max_ijk
    174     ENDIF
    175 
     227#endif
    176228
    177229    IF ( .NOT. dt_fixed )  THEN
     230#if defined( __openacc )
    178231!
    179232!--    Variable time step:
    180 !
    181 !--    For each component, compute the maximum time step according to the
    182 !--    CFL-criterion.
    183        dt_u = dx / ( ABS( u_max_cfl ) + 1.0E-10 )
    184        dt_v = dy / ( ABS( v_max_cfl ) + 1.0E-10 )
    185        dt_w = dzu(MAX( 1, w_max_ijk(1) )) / ( ABS( w_max ) + 1.0E-10 )
     233!--    Calculate the maximum time step according to the CFL-criterion,
     234!--    individually for each velocity component
     235       dt_u_l  =  999999.9
     236       dt_v_l  =  999999.9
     237       dt_w_l  =  999999.9
     238       u_max_l = -999999.9
     239       u_min_l =  999999.9
     240       v_max_l = -999999.9
     241       v_min_l =  999999.9
     242       w_max_l = -999999.9
     243       w_min_l =  999999.9
     244       !$acc parallel loop collapse(3) present( u, v, w )
     245       DO  i = nxl, nxr
     246          DO  j = nys, nyn
     247             DO  k = nzb+1, nzt
     248                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )
     249                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )
     250                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10 ) ) )
     251                u_max_l = MAX( u_max_l, u(k,j,i) )
     252                u_min_l = MIN( u_min_l, u(k,j,i) )
     253                v_max_l = MAX( v_max_l, v(k,j,i) )
     254                v_min_l = MIN( v_min_l, v(k,j,i) )
     255                w_max_l = MAX( w_max_l, w(k,j,i) )
     256                w_min_l = MIN( w_min_l, w(k,j,i) )
     257             ENDDO
     258          ENDDO
     259       ENDDO
     260       !$acc end parallel
     261
     262#if defined( __parallel )
     263       reduce_l(1) = dt_u_l
     264       reduce_l(2) = dt_v_l
     265       reduce_l(3) = dt_w_l
     266       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     267       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
     268       dt_u = reduce(1)
     269       dt_v = reduce(2)
     270       dt_w = reduce(3)
     271
     272       reduce_l(1) = u_max_l
     273       reduce_l(2) = v_max_l
     274       reduce_l(3) = w_max_l
     275       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     276       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
     277       u_max = reduce(1)
     278       v_max = reduce(2)
     279       w_max = reduce(3)
     280       reduce_l(1) = u_min_l
     281       reduce_l(2) = v_min_l
     282       reduce_l(3) = w_min_l
     283       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     284       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
     285       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
     286       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
     287       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
     288#else
     289       dt_u = dt_u_l
     290       dt_v = dt_v_l
     291       dt_w = dt_w_l
     292
     293       IF ( ABS( u_min_l ) > u_max_l )  THEN
     294          u_max = u_min_l
     295       ELSE
     296          u_max = u_max_l
     297       ENDIF
     298       IF ( ABS( v_min_l ) > v_max_l )  THEN
     299          v_max = v_min_l
     300       ELSE
     301          v_max = v_max_l
     302       ENDIF
     303       IF ( ABS( w_min_l ) > w_max_l )  THEN
     304          w_max = w_min_l
     305       ELSE
     306          w_max = w_max_l
     307       ENDIF
     308#endif
     309
     310#else
     311!
     312!--    Variable time step:
     313!--    Calculate the maximum time step according to the CFL-criterion,
     314!--    individually for each velocity component
     315       dt_u_l = 999999.9
     316       dt_v_l = 999999.9
     317       dt_w_l = 999999.9
     318       DO  i = nxl, nxr
     319          DO  j = nys, nyn
     320             DO  k = nzb+1, nzt
     321                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )
     322                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )
     323                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10 ) ) )
     324             ENDDO
     325          ENDDO
     326       ENDDO
     327
     328#if defined( __parallel )
     329       reduce_l(1) = dt_u_l
     330       reduce_l(2) = dt_v_l
     331       reduce_l(3) = dt_w_l
     332       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     333       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
     334       dt_u = reduce(1)
     335       dt_v = reduce(2)
     336       dt_w = reduce(3)
     337#else
     338       dt_u = dt_u_l
     339       dt_v = dt_v_l
     340       dt_w = dt_w_l
     341#endif
     342
     343#endif
    186344
    187345!
     
    199357!$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
    200358!$OMP DO
     359       !$acc parallel loop collapse(3) present( kh, km )
    201360       DO  i = nxl, nxr
    202361          DO  j = nys, nyn
    203362             DO  k = nzb+1, nzt
    204                 value = dxyz2_min(k) / ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 )
    205 
    206                 dt_diff_l = MIN( value, dt_diff_l )
     363                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / &
     364                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 ) )
    207365             ENDDO
    208366          ENDDO
    209367       ENDDO
     368       !$acc end parallel
    210369!$OMP END PARALLEL
    211370#if defined( __parallel )
     
    334493               '&dt_diff         = ', dt_diff, ' s',                          &
    335494               '&dt_plant_canopy = ', dt_plant_canopy, ' s',                  &
    336                '&u_max_cfl   = ', u_max_cfl, ' m/s   k=', u_max_cfl_ijk(1),   &
     495               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),       &
    337496               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
    338                '&v_max_cfl   = ', v_max_cfl, ' m/s   k=', v_max_cfl_ijk(1),   &
     497               '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),       &
    339498               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
    340                '&w_max       = ', w_max, ' m/s   k=', w_max_ijk(1),           &
     499               '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),       &
    341500               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
    342501          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
Note: See TracChangeset for help on using the changeset viewer.