subroutine nearwall
use modglobal, only:ib, ie, jb, je, jgb, jge, jmax, kb, ke, xh, xf, dy, zh, zf, lwarmstart, nblocks, libm, lzerogradtop, lwalldist
use modsubgriddata, only:lsmagorinsky, loneeqn
use modfields, only:mindist, wall
use modibmdata, only:xwallsglobal, ywallsglobal, zwallsglobal
use modmpi, only:myid
use initfac, only:block
implicit none
integer, allocatable :: ux0all(:, :, :), vy0all(:, :, :), wz0all(:, :, :)
real, allocatable :: distxf(:, :), distxh(:, :), distyf(:, :), distyh(:, :), distzf(:, :), distzh(:, :), &
distxf2(:, :), distxh2(:, :), distyf2(:, :), distyh2(:, :), distzf2(:, :), distzh2(:, :), distance(:)
real distx, disty, distz ! distx is the distance to nearest x-wall, etc.
! integer, allocatable :: optie(:)
integer ic, jc, kc, i, j, k, optie, il, iu, jl, ju, kl, ku, n
! if (lwarmstart .or. lles.eqv..false. .or. lvreman) then
if (((lsmagorinsky) .or. (loneeqn)) .and. (lwalldist)) then
if (myid == 0) then
write (6, *) 'Computing wall distances'
end if
! if (lles.eqv..false. .or. lvreman) then
mindist = 1.0e10
allocate (ux0all(ib - 1:ie + 1, jgb - 1:jge + 1, kb - 1:ke + 1)) ! This contains ux0 + the lower and (possibly) the upper wall
allocate (vy0all(ib - 1:ie + 1, jgb - 1:jge + 1, kb - 1:ke + 1)) ! This contains ux0 + the lower and (possibly) the upper wall
allocate (wz0all(ib - 1:ie + 1, jgb - 1:jge + 1, kb - 1:ke + 1)) ! This contains ux0 + the lower and (possibly) the upper wall
allocate (distxh(ib:ie, ib:ie + 1))
allocate (distxf(ib:ie, ib:ie + 1))
allocate (distyh(jb:je, jgb:jge + 1))
allocate (distyf(jb:je, jgb:jge + 1))
allocate (distzh(kb:ke, kb:ke + 1))
allocate (distzf(kb:ke, kb:ke + 1))
allocate (distxh2(ib:ie, ib:ie + 1))
allocate (distxf2(ib:ie, ib:ie + 1))
allocate (distyh2(jb:je, jgb:jge + 1))
allocate (distyf2(jb:je, jgb:jge + 1))
allocate (distzh2(kb:ke, kb:ke + 1))
allocate (distzf2(kb:ke, kb:ke + 1))
allocate (distance(4))
! initialize wall indicators
ux0all = 0
vy0all = 0
wz0all = 0
! Determine for each cell face if an x/y/z-wall is present
! from immersed boundaries
if (libm) then
! do loop over blocks
do n = 1, nblocks
il = block(n, 1)
iu = block(n, 2)
jl = block(n, 3)
ju = block(n, 4)
kl = block(n, 5)
ku = block(n, 6)
do k = kl, ku
do j = jl, ju
ux0all(il, j, k) = 1 ! lower x-wall
ux0all(iu + 1, j, k) = 1 ! upper x-wall
end do
end do
do k = kl, ku
do i = il, iu
vy0all(i, jl, k) = 1 ! lower y-wall
vy0all(i, ju + 1, k) = 1 ! upper y-wall
end do
end do
do j = jl, ju
do i = il, iu
wz0all(i, j, kl) = 1 ! lower z-wall
wz0all(i, j, ku + 1) = 1 ! upper z-wall
end do
end do
end do ! loop over nblocks
end if ! libm = .true.
! add the global walls (probably upper and lower wall, or only lower wall)
if (lzerogradtop) then
do i = ib, ie
do j = jgb, jge
wz0all(i, j, kb) = 1 ! ground wall
end do
end do
else
do i = ib, ie
do j = jgb, jge
wz0all(i, j, kb) = 1 ! ground wall
wz0all(i, j, ke + 1) = 1; ! top wall
end do
end do
end if
write (6, *) 'Determing distance matrices, proc=', myid
! Determine x-distance matrices:
do ic = ib, ie ! cell-center index
do i = ib, ie + 1 ! vertex-index (1 more than cell centers)
distxh(ic, i) = xf(ic) - xh(i)
end do
end do
do ic = ib, ie ! cell-center index
do i = ib, ie + 1 ! center -index
distxf(ic, i) = xf(ic) - xf(i)
end do
end do
! Determine y-distance matrices:
do jc = jb, je ! cell-center index
do j = jgb, jge + 1 ! vertex-index (1 more than cell centers) (global index to make sure distance to all cells is determined)
distyh(jc, j) = (jc + myid*jmax - j)*dy + 0.5*dy
end do
end do
do jc = jb, je ! cell-center index
do j = jgb, jge + 1 ! center-index (global index to make sure distance to all cells is determined)
distyf(jc, j) = (jc + myid*jmax - j)*dy
end do
end do
! Determine z-distance matrices:
do kc = kb, ke ! cell-center index
do k = kb, ke + 1 ! vertex-index (1 more than cell centers)
distzh(kc, k) = zf(kc) - zh(k)
end do
end do
do kc = kb, ke ! cell-center index
do k = kb, ke + 1 ! vertex-index (1 more than cell centers)
distzf(kc, k) = zf(kc) - zf(k)
end do
end do
distxh2 = distxh**2
distyh2 = distyh**2
distzh2 = distzh**2
distxf2 = distxf**2
distyf2 = distyf**2
distzf2 = distzf**2
write (6, *) 'Finished determing distance matrices, proc=', myid
write (6, *) 'determing distance to nearest wall for each cell center, proc=', myid
! Loop over cells (ic,jc,kc) for which minimal wall-to-cell-center-distance needs to be determined
! do jc=jgb,jge
do kc = kb, ke
do jc = jb, je
do ic = ib, ie
! Determine distance between cc of cell (ic,jc,kc) and faces of all cells (i,j,k)
do k = kb, ke + 1 ! Level ke+1 is computed in a separate loop (only necessary with upper wall-> global approach=faster)
do j = jgb, jge + 1 ! loop goes up to jge+1 because jge+1 contains the last vy0-wall
do i = ib, ie + 1 ! loop goes up to ie+1 because ie+1 contains the last ux0-wall
if (ux0all(i, j, k) == 1 .OR. vy0all(i, j, k) == 1 .OR. wz0all(i, j, k) == 1) then
distx = 1.0e10 ! make sure distx is very large when no x-wall is present
disty = 1.0e10 ! make sure disty is very large when no y-wall is present
distz = 1.0e10 ! make sure distz is very large when no z-wall is present
if (ux0all(i, j, k) == 1) then
distx = sqrt(distxh2(ic, i) + distyf2(jc, j) + distzf2(kc, k))
end if
if (vy0all(i, j, k) == 1) then
disty = sqrt(distxf2(ic, i) + distyh2(jc, j) + distzf2(kc, k))
end if
if (wz0all(i, j, k) == 1) then
distz = sqrt(distxf2(ic, i) + distyf2(jc, j) + distzh2(kc, k))
end if
else ! no walls are present in cell (i,j,k) -> distance does not need to be determined for this cell
cycle ! go to next cell (i,j,k)
end if
! determine minimal wall distance between cc of (ic,jc,kc) and faces of cell (i,j,k)
distance = (/mindist(ic, jc, kc), distx, disty, distz/)
optie = minloc(distance, 1)
! write(6,*) 'optie=', optie
if (optie == 1) then
cycle
else if (optie == 2) then
mindist(ic, jc, kc) = distx
wall(ic, jc, kc, 2) = j
wall(ic, jc, kc, 3) = k
! wall(ic,jc,kc,4) = 1 ! This means the wall closest to the cc of (ic,jc,kc) is at an x-wall at (i,j,k)
if (ic >= i) then
wall(ic, jc, kc, 1) = i
wall(ic, jc, kc, 4) = 5 ! shear component index (variable: shear)
wall(ic, jc, kc, 5) = 9 ! shear component index (variable: shear)
else
wall(ic, jc, kc, 1) = i - 1 ! in the subgrid this stress is computed in the cell i-1
wall(ic, jc, kc, 4) = 6 ! shear component index (variable: shear)
wall(ic, jc, kc, 5) = 10 ! shear component index (variable: shear)
end if
else if (optie == 3) then
mindist(ic, jc, kc) = disty
wall(ic, jc, kc, 1) = i
wall(ic, jc, kc, 3) = k
! wall(ic,jc,kc,4) = 2 ! This means the wall closest to the cc of (ic,jc,kc) is at a y-wall at (i,j,k)
if (jc + myid*jmax >= j) then
wall(ic, jc, kc, 2) = j
wall(ic, jc, kc, 4) = 1 ! shear component index (variable: shear)
wall(ic, jc, kc, 5) = 11 ! shear component index (variable: shear)
else
wall(ic, jc, kc, 2) = j - 1 ! in the subgrid this stress is computed in the cell j-1
wall(ic, jc, kc, 4) = 2 ! shear component index (variable: shear)
wall(ic, jc, kc, 5) = 12 ! shear component index (variable: shear)
end if
else if (optie == 4) then
mindist(ic, jc, kc) = distz
wall(ic, jc, kc, 1) = i
wall(ic, jc, kc, 2) = j
! wall(ic,jc,kc,4) = 3 ! This means the wall closest to the cc of (ic,jc,kc) is at a z-wall at (i,j,k)
if (kc >= k) then
wall(ic, jc, kc, 3) = k
wall(ic, jc, kc, 4) = 3 ! shear component index (variable: shear)
wall(ic, jc, kc, 5) = 7 ! shear component index (variable: shear)
else
wall(ic, jc, kc, 3) = k - 1 ! in the subgrid this stress is computed in the cel k-1
wall(ic, jc, kc, 4) = 4 ! shear component index (variable: shear)
wall(ic, jc, kc, 5) = 8 ! shear component index (variable: shear)
end if
end if
! mindist(ic,jc+myid*jmax,kc)=min(mindist(ic,jc+myid*jmax,kc),distx,disty,distz) ! global j index
end do
end do
end do
! if (myid==0) write(6,*) 'finished for cell (ic,jc,kc)=',ic,jc,kc
end do
end do
end do
write (6, *) 'Finished determing distance to nearest wall for each cell center, proc=', myid
! write(6,*) 'mindist(ib,jb+myid*jmax,kb),mindist(ib,je+myid*jmax,kb)',mindist(ib,jb+myid*jmax,kb),mindist(ib,je+myid*jmax,kb)
else
return
end if !(lwarmstart)
deallocate (ux0all, vy0all, wz0all)
deallocate (xwallsglobal, ywallsglobal, zwallsglobal, block) ! used for determining boundaries
return
end subroutine nearwall