Changeset 4101 for palm/trunk/SOURCE/timestep.f90

Timestamp:

Jul 17, 2019 3:14:26 PM (5 years ago)

Author:

gronemeier

Message:

timestep.f90:

• consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
• in RANS mode, instead of considering each wind component individually use the wind speed of 3d wind vector in CFL criterion
• do not limit the increase of dt based on its previous value in RANS mode

other:

• remove dt_old

File:

• palm/trunk/SOURCE/timestep.f90 (modified) (7 diffs)

palm/trunk/SOURCE/timestep.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! - consider 2*Km within diffusion criterion as Km is considered twice within
+!   the diffusion of e,
+! - in RANS mode, instead of considering each wind component individually use
+!   the wind speed of 3d wind vector in CFL criterion
+! - do not limit the increase of dt based on its previous value in RANS mode
+!
+! 3658 2019-01-07 20:28:54Z knoop
 ! OpenACC port for SPEC
 ! 
@@ -146 +153 @@
     USE control_parameters,                                                    &
         ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
-               galilei_transformation, old_dt, message_string, rans_mode,      &
+               galilei_transformation, message_string, rans_mode,              &
                stop_dt, terminate_coupled, terminate_coupled_remote,           &
                timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans
@@ -218 +225 @@
     !$ACC HOST(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
     !$ACC HOST(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
-
-!
-!--    Save former time step as reference
-       old_dt = dt_3d
 
 !
@@ -284 +287 @@
 !
 !--    Variable time step:
-!--    Calculate the maximum time step according to the CFL-criterion,
-!--    individually for each velocity component
+!--    Calculate the maximum time step according to the CFL-criterion
        dt_u_l = 999999.9_wp
        dt_v_l = 999999.9_wp
        dt_w_l = 999999.9_wp
-       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
-       !$ACC COPY(dt_u_l, dt_v_l, dt_w_l, u_stokes_zu, v_stokes_zu) &
-       !$ACC REDUCTION(MIN: dt_u_l, dt_v_l, dt_w_l) &
-       !$ACC PRESENT(u, v, w, dzu)
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                dt_u_l = MIN( dt_u_l, ( dx     /                               &
-                                 ( ABS( u(k,j,i) - u_gtrans + u_stokes_zu(k) ) &
-                                   + 1.0E-10_wp ) ) )
-                dt_v_l = MIN( dt_v_l, ( dy     /                               &
-                                 ( ABS( v(k,j,i) - v_gtrans + v_stokes_zu(k) ) &
-                                   + 1.0E-10_wp ) ) )
-                dt_w_l = MIN( dt_w_l, ( dzu(k) /                               &
-                                 ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
+
+       IF ( .NOT. rans_mode )  THEN
+!
+!--       Consider each velocity component individually
+
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC COPY(dt_u_l, dt_v_l, dt_w_l, u_stokes_zu, v_stokes_zu) &
+          !$ACC REDUCTION(MIN: dt_u_l, dt_v_l, dt_w_l) &
+          !$ACC PRESENT(u, v, w, dzu)
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   dt_u_l = MIN( dt_u_l, ( dx     /                               &
+                                    ( ABS( u(k,j,i) - u_gtrans + u_stokes_zu(k) ) &
+                                      + 1.0E-10_wp ) ) )
+                   dt_v_l = MIN( dt_v_l, ( dy     /                               &
+                                    ( ABS( v(k,j,i) - v_gtrans + v_stokes_zu(k) ) &
+                                      + 1.0E-10_wp ) ) )
+                   dt_w_l = MIN( dt_w_l, ( dzu(k) /                               &
+                                    ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
+                ENDDO
              ENDDO
           ENDDO
-       ENDDO
+
+       ELSE
+!
+!--       Consider the wind speed at the scalar-grid point
+!--       !> @note considering the wind speed instead of each individual wind
+!--       !>       component is only a workaround so far. This might has to be
+!--       !>       changed in the future.
+
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC COPY(dt_u_l, u_stokes_zu, v_stokes_zu) &
+          !$ACC REDUCTION(MIN: dt_u_l) &
+          !$ACC PRESENT(u, v, w, dzu)
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   dt_u_l = MIN( dt_u_l, ( MIN( dx, dy, dzu(k) ) / ( &
+                      SQRT(  ( 0.5 * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans + u_stokes_zu(k) )**2   &
+                           + ( 0.5 * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans + v_stokes_zu(k) )**2   &
+                           + ( 0.5 * ( w(k,j,i) + w(k-1,j,i) )                             )**2 ) &
+                      + 1.0E-10_wp ) ) )
+                ENDDO
+             ENDDO
+          ENDDO
+          
+          dt_v_l = dt_u_l
+          dt_w_l = dt_u_l
+
+       ENDIF
 
 #if defined( __parallel )
@@ -342 +377 @@
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
-                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) /                     &
-                                    ( MAX( kh(k,j,i), ABS( km(k,j,i) ) )       &
+                dt_diff_l = MIN( dt_diff_l,                                       &
+                                 dxyz2_min(k) /                                   &
+                                    ( MAX( kh(k,j,i), 2.0_wp * ABS( km(k,j,i) ) ) &
                                       + 1E-20_wp ) )
              ENDDO
@@ -363 +399 @@
        dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
        dt_3d = MIN( dt_3d, dt_max )
-!
-!--    In RANS mode, the time step must not increase by more than a factor of 2
-       IF ( rans_mode )  dt_3d = MIN( dt_3d, dt_max, 2.0_wp * old_dt )
 
 !
@@ -390 +423 @@
           WRITE( message_string, * ) 'Time step has reached minimum limit.',   &
                '&dt              = ', dt_3d, ' s  Simulation is terminated.',  &
-               '&old_dt          = ', old_dt, ' s',                            &
                '&dt_u            = ', dt_u, ' s',                              &
                '&dt_v            = ', dt_v, ' s',                              &