fix computation of retau for pdc case (CaNS-World#64)

* fix computation of retau for `pdc` case (it only would work for a case with half channel height `lz/2 = 1`)

* added a constant pressure gradient turbulent half-channel case (`hdc`).

* avoid `sqrt` of negative values for means near zero under `output.f90` by computing the variance rather than the r.m.s.

examples/turbulent_half_channel_constant_pressure_gradient/dns.in

@@ -0,0 +1,23 @@
+512 256 72               ! itot, jtot, ktot
+12. 6. 1.                ! lx, ly, lz
+0.                       ! gr
+.95 1.e5                 ! cfl, dtmin
+180.                     ! visci
+hdc                      ! inivel
+T                        ! is_wallturb
+100000 100. 0.1          ! nstep,time_max,tw_max
+T F F                    ! stop_type(1:3)
+F T 0                    ! restart,is_overwrite_save,nsaves_max
+10 10 100 500 10000 5000 ! icheck, iout0d, iout1d, iout2d, iout3d, isave
+P P  P P  D D            ! cbcvel(0:1,1:3,1) [u BC type]
+P P  P P  D D            ! cbcvel(0:1,1:3,2) [v BC type]
+P P  P P  D D            ! cbcvel(0:1,1:3,3) [w BC type]
+P P  P P  N N            ! cbcpre(0:1,1:3  ) [p BC type]
+0. 0.  0. 0.  0. 0.      !  bcvel(0:1,1:3,1) [u BC value]
+0. 0.  0. 0.  0. 0.      !  bcvel(0:1,1:3,2) [v BC value]
+0. 0.  0. 0.  0. 0.      !  bcvel(0:1,1:3,3) [w BC value]
+0. 0.  0. 0.  0. 0.      !  bcpre(0:1,1:3  ) [p BC value]
+1. 0. 0.                 ! bforce(1:3)
+F F F                    ! is_forced(1:3)
+0. 0. 0.                 ! velf(1:3)
+2 2                      ! dims(1:2)

src/initflow.f90

@@ -29,7 +29,7 @@ subroutine initflow(inivel,ng,lo,zc,dzc,dzf,visc,u,v,w,p)
     logical :: is_noise,is_mean,is_pair
     real(rp) :: xc,yc,zcc,xf,yf,zff
     real(rp), allocatable, dimension(:) :: zc2
-    real(rp) :: uref
+    real(rp) :: uref,lref
     real(rp) :: ubulk,reb,retau
     integer, dimension(3) :: n
     !
@@ -121,16 +121,29 @@ subroutine initflow(inivel,ng,lo,zc,dzc,dzf,visc,u,v,w,p)
           end do
         end do
       end do
-    case('pdc')
+    case('pdc','hdc')
+      lref  = lz/2.
+      if(trim(inivel) /= 'pdc') lref = 2.*lref
       if(is_wallturb) then ! turbulent flow
-        uref  = (bforce(1)/(lz/2.))**(0.5)
-        retau = uref*(lz/2.)/visc
+        uref  = (bforce(1)*lref)**(0.5) ! utau = sqrt(-dpdx*h)
+        retau = uref*lref/visc
         reb   = (retau/.09)**(1./.88)
-        ubulk = (reb/2.)/retau*uref
+        ubulk = reb*visc/(2*lref)
       else                 ! laminar flow
-        ubulk = (bforce(1)*(lz/2.)**2/(3.*visc))
+        ubulk = (bforce(1)*lref**2/(3.*visc))
+      end if
+      if(trim(inivel) == 'pdc') then
+        call poiseuille(n(3),zc/lz,ubulk,u1d)
+      else
+        deallocate(u1d)
+        allocate(u1d(2*n(3)))
+        allocate(zc2(0:2*n(3)+1))
+        zc2(1     :  n(3)) =        zc(1   :n(3): 1)
+        zc2(n(3)+1:2*n(3)) = 2*lz - zc(n(3):1   :-1)
+        zc2(0)        = -zc(0)
+        zc2(2*n(3)+1) = 2*lz + zc(0)
+        call poiseuille(2*n(3),zc2/(2*lz),ubulk,u1d)
       end if
-      call poiseuille(n(3),zc/lz,ubulk,u1d)
       is_mean=.true.
     case default
       if(myid == 0) print*, 'ERROR: invalid name for initial velocity field'

src/output.f90

@@ -349,9 +349,9 @@ subroutine out1d_chan(fname,ng,lo,hi,idir,l,dl,z_g,u,v,w) ! e.g. for a channel w
       um(:) = um(:)*grid_area_ratio
       vm(:) = vm(:)*grid_area_ratio
       wm(:) = wm(:)*grid_area_ratio
-      u2(:) = sqrt(u2(:)*grid_area_ratio - um(:)**2)
-      v2(:) = sqrt(v2(:)*grid_area_ratio - vm(:)**2)
-      w2(:) = sqrt(w2(:)*grid_area_ratio - wm(:)**2)
+      u2(:) = u2(:)*grid_area_ratio - um(:)**2
+      v2(:) = v2(:)*grid_area_ratio - vm(:)**2
+      w2(:) = w2(:)*grid_area_ratio - wm(:)**2
       uw(:) = uw(:)*grid_area_ratio - um(:)*wm(:)
       if(myid == 0) then
         open(newunit=iunit,file=fname)