subroutine initibmwallfun(fname_bnd, fname_sec, dir, bound_info)
use modglobal, only : ifinput, ib, itot, ih, jb, jtot, jh, kb, ke, ktot, kh, &
xf, yf, zf, xh, yh, zh, dx, dy, dzf, xhat, yhat, zhat, eps1
use modmpi, only : myid, comm3d, MY_REAL, mpierr
use initfac, only : facnorm, facz0
use decomp_2d, only : zstart, zend
use readinput, only : read_sparse_ijk
character(20), intent(in) :: fname_bnd, fname_sec
type(bound_info_type) :: bound_info
real, intent(in), dimension(3) :: dir
real, dimension(ib:itot+ih) :: xgrid
real, dimension(jb:jtot+jh) :: ygrid
real, dimension(kb:ktot+kh) :: zgrid
logical, dimension(bound_info%nfctsecs) :: lfctsecsrank
integer, parameter :: nchunk = 262144 ! sections per broadcast chunk (~20 MB of buffers)
integer :: c0, nc, nn
integer, allocatable :: ibuf1(:), ibuf2(:), ibufu(:,:), ibufv(:,:), ibufw(:,:), ibufc(:,:)
real, allocatable :: rbuf1(:), rbuf2(:), rbuf3(:,:)
logical, allocatable :: lbuf1(:), lbuf2(:)
logical, dimension(:), allocatable :: lbndptsrank
real, dimension(3) :: norm, p0, p1, pxl, pxu, pyl, pyu, pzl, pzu
integer, dimension(6) :: check
integer, dimension(1) :: pos_min_dist
real, dimension(6,3) :: inter
real, dimension(6) :: inter_dists
real :: xc, yc, zc, xl, yl, zl, xu, yu, zu
integer n, m, norm_align, dir_align, pos, p
character(80) :: chmess
integer, dimension(:), allocatable :: ids_loc
integer, dimension(:,:), allocatable :: pts_loc
! Read boundary points using generic read_sparse_ijk routine (skips 1 header line)
! Request both local and global arrays since sections need global indices
call read_sparse_ijk(fname_bnd, bound_info%nbndpts, bound_info%nbndptsrank, ids_loc, pts_loc, nskip=1, pts_glob_out=bound_info%bndpts)
! Transfer ownership of local arrays (no copying, no conversion needed)
call move_alloc(ids_loc, bound_info%bndptsrank)
call move_alloc(pts_loc, bound_info%bndpts_loc)
! Build lbndptsrank lookup array for determining which sections are on this rank
allocate(lbndptsrank(bound_info%nbndpts))
lbndptsrank = .false.
allocate(ibuf1(nchunk), ibuf2(nchunk), rbuf1(nchunk), rbuf2(nchunk), rbuf3(nchunk,3), &
ibufu(nchunk,3), ibufv(nchunk,3), ibufw(nchunk,3), ibufc(nchunk,3), lbuf1(nchunk), lbuf2(nchunk))
do m = 1, bound_info%nbndptsrank
lbndptsrank(bound_info%bndptsrank(m)) = .true.
end do
! The global (all-sections) arrays are needed only on rank 0, which reads the file
! and does the reconstruction geometry; they are then handed out in chunks so that
! the other ranks never hold nfctsecs-sized arrays (about 1 GB per rank for a
! city-scale case, and the start-up memory peak of the model).
if (myid == 0) then
allocate(bound_info%secfacids(bound_info%nfctsecs))
allocate(bound_info%secareas(bound_info%nfctsecs))
allocate(bound_info%secbndptids(bound_info%nfctsecs))
allocate(bound_info%bnddst(bound_info%nfctsecs))
allocate(bound_info%recpts(bound_info%nfctsecs,3))
allocate(bound_info%recids_u(bound_info%nfctsecs,3))
allocate(bound_info%recids_v(bound_info%nfctsecs,3))
allocate(bound_info%recids_w(bound_info%nfctsecs,3))
allocate(bound_info%recids_c(bound_info%nfctsecs,3))
allocate(bound_info%lcomprec(bound_info%nfctsecs))
allocate(bound_info%lskipsec(bound_info%nfctsecs))
end if
dir_align = alignment(dir)
select case(dir_align)
case(1)
xgrid = xh
ygrid = yf
zgrid = zf
case(2)
xgrid = xf
ygrid = yh
zgrid = zf
case(3)
xgrid = xf
ygrid = yf
zgrid = zh
case(0)
xgrid = xf
ygrid = yf
zgrid = zf
end select
if (myid == 0) then
open (ifinput, file=fname_sec)
read (ifinput, '(a80)') chmess
do n = 1, bound_info%nfctsecs
read (ifinput, *) bound_info%secfacids(n), bound_info%secareas(n), bound_info%secbndptids(n), bound_info%bnddst(n)
!bound_info%intpts(n,1), bound_info%intpts(n,2), bound_info%intpts(n,3)
end do
close (ifinput)
do n = 1,bound_info%nfctsecs
m = bound_info%secbndptids(n)
!bound_info%bndvec(n,1) = xgrid(bound_info%bndpts(m,1)) - bound_info%intpts(n,1)
!bound_info%bndvec(n,2) = ygrid(bound_info%bndpts(m,2)) - bound_info%intpts(n,2)
!bound_info%bndvec(n,3) = zgrid(bound_info%bndpts(m,3)) - bound_info%intpts(n,3)
!bound_info%bnddst(n) = norm2(bound_info%bndvec(n,:))
!write(*,*) bound_info%bnddst(n)
!bound_info%bndvec(n,:) = bound_info%bndvec(n,:) / bound_info%bnddst(n)
norm = facnorm(bound_info%secfacids(n),:)
norm_align = alignment(norm)
if ((dir_align /= 0 .and. dir_align == norm_align) .or. (facz0(bound_info%secfacids(n)) < eps1)) then
! (for velocities) if the facet is aligned with the grid AND in the same direction as the current velocity grid direction
! therefore no tangential component, don't need to calculate shear stress
bound_info%lskipsec(n) = .true.
cycle
else
bound_info%lskipsec(n) = .false.
end if
if (log(bound_info%bnddst(n)/facz0(bound_info%secfacids(n))) > 1. .or. lnorec) then ! the wall function is well-defined
bound_info%lcomprec(n) = .true. ! do simple reconstruction
else ! need to reconstruct
bound_info%lcomprec(n) = .false.
! Find reconstruction point
! cell centre (of current grid)
xc = xgrid(bound_info%bndpts(m,1))
yc = ygrid(bound_info%bndpts(m,2))
zc = zgrid(bound_info%bndpts(m,3))
! cell edges
xl = xc - dx/2.
xu = xc + dx/2.
yl = yc - dy/2.
yu = yc + dy/2.
zl = zc - dzf(1)/2. ! assumes equidistant
zu = zc + dzf(1)/2. ! assumes equidistant
! points on planes
pxl = (/xl, yc, zc/)
pxu = (/xu, yc, zc/)
pyl = (/xc, yl, zc/)
pyu = (/xc, yu, zc/)
pzl = (/xc, yc, zl/)
pzu = (/xc, yc, zu/)
p0 = (/xc, yc, zc/)
p1 = p0 + norm * sqrt(3.)*(dx*dy*dzf(1))**(1./3.)
call plane_line_intersection(xhat, pxl, p0, p1, inter(1,:), check(1), inter_dists(1))
call plane_line_intersection(xhat, pxu, p0, p1, inter(2,:), check(2), inter_dists(2))
call plane_line_intersection(yhat, pyl, p0, p1, inter(3,:), check(3), inter_dists(3))
call plane_line_intersection(yhat, pyu, p0, p1, inter(4,:), check(4), inter_dists(4))
call plane_line_intersection(zhat, pzl, p0, p1, inter(5,:), check(5), inter_dists(5))
call plane_line_intersection(zhat, pzu, p0, p1, inter(6,:), check(6), inter_dists(6))
pos_min_dist = minloc(inter_dists, mask=check==1)
pos = pos_min_dist(1)
if (pos == 0) then
write(*,*) "ERROR: no intersection found"
stop 1
else
bound_info%recpts(n,:) = inter(pos,:) ! x y z
end if
! find which cell the point lies in
bound_info%recids_u(n,1) = findloc(bound_info%recpts(n,1) >= xh, .true., 1, back=.true.)
bound_info%recids_u(n,2) = findloc(bound_info%recpts(n,2) >= yf, .true., 1, back=.true.)
bound_info%recids_u(n,3) = findloc(bound_info%recpts(n,3) >= zf, .true., 1, back=.true.)
bound_info%recids_v(n,1) = findloc(bound_info%recpts(n,1) >= xf, .true., 1, back=.true.)
bound_info%recids_v(n,2) = findloc(bound_info%recpts(n,2) >= yh, .true., 1, back=.true.)
bound_info%recids_v(n,3) = findloc(bound_info%recpts(n,3) >= zf, .true., 1, back=.true.)
bound_info%recids_w(n,1) = findloc(bound_info%recpts(n,1) >= xf, .true., 1, back=.true.)
bound_info%recids_w(n,2) = findloc(bound_info%recpts(n,2) >= yf, .true., 1, back=.true.)
bound_info%recids_w(n,3) = findloc(bound_info%recpts(n,3) >= zh, .true., 1, back=.true.)
bound_info%recids_c(n,1) = findloc(bound_info%recpts(n,1) >= xf, .true., 1, back=.true.)
bound_info%recids_c(n,2) = findloc(bound_info%recpts(n,2) >= yf, .true., 1, back=.true.)
bound_info%recids_c(n,3) = findloc(bound_info%recpts(n,3) >= zf, .true., 1, back=.true.)
! check to see if recids is inside the domain
if (bound_info%recids_u(n,1) < ib .or. bound_info%recids_u(n,1)+1 > itot+ih .or. &
bound_info%recids_u(n,2) < jb .or. bound_info%recids_u(n,2)+1 > jtot+jh .or. &
bound_info%recids_u(n,3) < kb .or. bound_info%recids_u(n,3)+1 > ke+kh) then
! if (myid == 0) then
! write(*,*) "DEBUG: skipping section n=", n, "- u-velocity reconstruction index out of bounds:", bound_info%recids_u(n,1), bound_info%recids_u(n,1)+1, bound_info%recids_u(n,2), bound_info%recids_u(n,2)+1, bound_info%recids_u(n,3), bound_info%recids_u(n,3)+1
! end if
bound_info%lskipsec(n) = .true.
cycle
end if
if (bound_info%recids_v(n,1) < ib .or. bound_info%recids_v(n,1)+1 > itot+ih .or. &
bound_info%recids_v(n,2) < jb .or. bound_info%recids_v(n,2)+1 > jtot+jh .or. &
bound_info%recids_v(n,3) < kb .or. bound_info%recids_v(n,3)+1 > ke+kh) then
! if (myid == 0) then
! write(*,*) "DEBUG: skipping section n=", n, "- v-velocity reconstruction index out of bounds:", bound_info%recids_v(n,1), bound_info%recids_v(n,1)+1, bound_info%recids_v(n,2), bound_info%recids_v(n,2)+1, bound_info%recids_v(n,3), bound_info%recids_v(n,3)+1
! end if
bound_info%lskipsec(n) = .true.
cycle
end if
if (bound_info%recids_w(n,1) < ib .or. bound_info%recids_w(n,1)+1 > itot+ih .or. &
bound_info%recids_w(n,2) < jb .or. bound_info%recids_w(n,2)+1 > jtot+jh .or. &
bound_info%recids_w(n,3) < kb .or. bound_info%recids_w(n,3)+1 > ke+kh) then
! if (myid == 0) then
! write(*,*) "DEBUG: skipping section n=", n, "- w-velocity reconstruction index out of bounds:", bound_info%recids_w(n,1), bound_info%recids_w(n,1)+1, bound_info%recids_w(n,2), bound_info%recids_w(n,2)+1, bound_info%recids_w(n,3), bound_info%recids_w(n,3)+1
! end if
bound_info%lskipsec(n) = .true.
cycle
end if
if (bound_info%recids_c(n,1) < ib .or. bound_info%recids_c(n,1)+1 > itot+ih .or. &
bound_info%recids_c(n,2) < jb .or. bound_info%recids_c(n,2)+1 > jtot+jh .or. &
bound_info%recids_c(n,3) < kb .or. bound_info%recids_c(n,3)+1 > ke+kh) then
! if (myid == 0) then
! write(*,*) "DEBUG: skipping section n=", n, "- c-velocity index out of bounds:", bound_info%recids_c(n,1), bound_info%recids_c(n,1)+1, bound_info%recids_c(n,2), bound_info%recids_c(n,2)+1, bound_info%recids_c(n,3), bound_info%recids_c(n,3)+1
! end if
bound_info%lskipsec(n) = .true.
cycle
end if
!check recpts is inside the box defined by the corners
! u
if ((bound_info%recpts(n,1) < xh(bound_info%recids_u(n,1))) .or. &
(bound_info%recpts(n,1) > xh(bound_info%recids_u(n,1)+1))) then
write(*,*) "ERROR: x out of bounds"
stop 1
end if
if ((bound_info%recpts(n,2) < yf(bound_info%recids_u(n,2))) .or. &
(bound_info%recpts(n,2) > yf(bound_info%recids_u(n,2)+1))) then
write(*,*) "ERROR: y out of bounds"
stop 1
end if
if ((bound_info%recpts(n,3) < zf(bound_info%recids_u(n,3))) .or. &
(bound_info%recpts(n,3) > zf(bound_info%recids_u(n,3)+1))) then
write(*,*) "ERROR: z out of bounds"
stop 1
end if
! v
if ((bound_info%recpts(n,1) < xf(bound_info%recids_v(n,1))) .or. &
(bound_info%recpts(n,1) > xf(bound_info%recids_v(n,1)+1))) then
write(*,*) "ERROR: x out of bounds"
stop 1
end if
if ((bound_info%recpts(n,2) < yh(bound_info%recids_v(n,2))) .or. &
(bound_info%recpts(n,2) > yh(bound_info%recids_v(n,2)+1))) then
write(*,*) "ERROR: y out of bounds"
stop 1
end if
if ((bound_info%recpts(n,3) < zf(bound_info%recids_v(n,3))) .or. &
(bound_info%recpts(n,3) > zf(bound_info%recids_v(n,3)+1))) then
write(*,*) "ERROR: z out of bounds"
stop 1
end if
! w
if ((bound_info%recpts(n,1) < xf(bound_info%recids_w(n,1))) .or. &
(bound_info%recpts(n,1) > xf(bound_info%recids_w(n,1)+1))) then
write(*,*) "ERROR: x out of bounds"
stop 1
end if
if ((bound_info%recpts(n,2) < yf(bound_info%recids_w(n,2))) .or. &
(bound_info%recpts(n,2) > yf(bound_info%recids_w(n,2)+1))) then
write(*,*) "ERROR: y out of bounds"
stop 1
end if
if ((bound_info%recpts(n,3) < zh(bound_info%recids_w(n,3))) .or. &
(bound_info%recpts(n,3) > zh(bound_info%recids_w(n,3)+1))) then
write(*,*) "ERROR: z out of bounds"
stop 1
end if
end if
end do
end if ! myid==0
! Pass 1: hand out the section -> boundary-point ids in chunks; a section belongs to
! this rank when its boundary point does. (lfctsecsrank is a logical per section.)
bound_info%nfctsecsrank = 0
do c0 = 1, bound_info%nfctsecs, nchunk
nc = min(nchunk, bound_info%nfctsecs - c0 + 1)
if (myid == 0) ibuf1(1:nc) = bound_info%secbndptids(c0:c0+nc-1)
call MPI_BCAST(ibuf1, nc, MPI_INTEGER, 0, comm3d, mpierr)
do n = 1, nc
if (lbndptsrank(ibuf1(n))) then
lfctsecsrank(c0+n-1) = .true.
bound_info%nfctsecsrank = bound_info%nfctsecsrank + 1
else
lfctsecsrank(c0+n-1) = .false.
end if
end do
end do
allocate(bound_info%fctsecsrank(bound_info%nfctsecsrank))
allocate(bound_info%secfacids_loc(bound_info%nfctsecsrank))
allocate(bound_info%secareas_loc(bound_info%nfctsecsrank))
allocate(bound_info%secbndpts_loc(bound_info%nfctsecsrank,3))
allocate(bound_info%bnddst_loc(bound_info%nfctsecsrank))
allocate(bound_info%recpts_loc(bound_info%nfctsecsrank,3))
allocate(bound_info%recids_u_loc(bound_info%nfctsecsrank,3))
allocate(bound_info%recids_v_loc(bound_info%nfctsecsrank,3))
allocate(bound_info%recids_w_loc(bound_info%nfctsecsrank,3))
allocate(bound_info%recids_c_loc(bound_info%nfctsecsrank,3))
allocate(bound_info%lcomprec_loc(bound_info%nfctsecsrank))
allocate(bound_info%lskipsec_loc(bound_info%nfctsecsrank))
! Pass 2: hand out all section data in the same chunks; each rank keeps its own rows
! in ascending section order, exactly as the former full-array selection did.
m = 0
p = 0 ! counter for sections overridden to simple reconstruction
do c0 = 1, bound_info%nfctsecs, nchunk
nc = min(nchunk, bound_info%nfctsecs - c0 + 1)
if (myid == 0) then
ibuf1(1:nc) = bound_info%secfacids(c0:c0+nc-1)
ibuf2(1:nc) = bound_info%secbndptids(c0:c0+nc-1)
rbuf1(1:nc) = bound_info%secareas(c0:c0+nc-1)
rbuf2(1:nc) = bound_info%bnddst(c0:c0+nc-1)
rbuf3(1:nc,:) = bound_info%recpts(c0:c0+nc-1,:)
ibufu(1:nc,:) = bound_info%recids_u(c0:c0+nc-1,:)
ibufv(1:nc,:) = bound_info%recids_v(c0:c0+nc-1,:)
ibufw(1:nc,:) = bound_info%recids_w(c0:c0+nc-1,:)
ibufc(1:nc,:) = bound_info%recids_c(c0:c0+nc-1,:)
lbuf1(1:nc) = bound_info%lcomprec(c0:c0+nc-1)
lbuf2(1:nc) = bound_info%lskipsec(c0:c0+nc-1)
end if
call MPI_BCAST(ibuf1, nc, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(ibuf2, nc, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(rbuf1, nc, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(rbuf2, nc, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(rbuf3, size(rbuf3), MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(ibufu, size(ibufu), MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(ibufv, size(ibufv), MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(ibufw, size(ibufw), MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(ibufc, size(ibufc), MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(lbuf1, nc, MPI_LOGICAL, 0, comm3d, mpierr)
call MPI_BCAST(lbuf2, nc, MPI_LOGICAL, 0, comm3d, mpierr)
do nn = 1, nc
n = c0 + nn - 1
if (.not. lfctsecsrank(n)) cycle
m = m + 1
bound_info%fctsecsrank(m) = n
bound_info%secfacids_loc(m) = ibuf1(nn) ! facet id
bound_info%secareas_loc(m) = rbuf1(nn)
bound_info%secbndpts_loc(m,:) = bound_info%bndpts(ibuf2(nn),:) ! boundary point location (in global coordinates)
if (bound_info%bndpts(ibuf2(nn),1) < zstart(1) .or. bound_info%bndpts(ibuf2(nn),1) > zend(1)) then
write(*,*) "problem in x boundary points on : ", myid, n, ibuf2(nn), bound_info%bndpts(ibuf2(nn),1), zstart(1), zend(1)
end if
if (bound_info%bndpts(ibuf2(nn),2) < zstart(2) .or. bound_info%bndpts(ibuf2(nn),2) > zend(2)) then
write(*,*) "problem in y boundary points on rank: ", myid, n, ibuf2(nn), bound_info%bndpts(ibuf2(nn),2), zstart(2), zend(2)
end if
bound_info%bnddst_loc(m) = rbuf2(nn)
bound_info%recpts_loc(m,:) = rbuf3(nn,:)
bound_info%recids_u_loc(m,:) = ibufu(nn,:)
bound_info%recids_v_loc(m,:) = ibufv(nn,:)
bound_info%recids_w_loc(m,:) = ibufw(nn,:)
bound_info%recids_c_loc(m,:) = ibufc(nn,:)
bound_info%lcomprec_loc(m) = lbuf1(nn)
if (.not. bound_info%lcomprec_loc(m) .and. .not. lbuf2(nn)) then
if (bound_info%recids_u_loc(m,1) < zstart(1)-1 .or. bound_info%recids_u_loc(m,1) > zend(1)+1 .or. &
bound_info%recids_u_loc(m,2) < zstart(2)-1 .or. bound_info%recids_u_loc(m,2) > zend(2)+1 .or. &
bound_info%recids_v_loc(m,1) < zstart(1)-1 .or. bound_info%recids_v_loc(m,1) > zend(1)+1 .or. &
bound_info%recids_v_loc(m,2) < zstart(2)-1 .or. bound_info%recids_v_loc(m,2) > zend(2)+1 .or. &
bound_info%recids_w_loc(m,1) < zstart(1)-1 .or. bound_info%recids_w_loc(m,1) > zend(1)+1 .or. &
bound_info%recids_w_loc(m,2) < zstart(2)-1 .or. bound_info%recids_w_loc(m,2) > zend(2)+1 .or. &
bound_info%recids_c_loc(m,1) < zstart(1)-1 .or. bound_info%recids_c_loc(m,1) > zend(1)+1 .or. &
bound_info%recids_c_loc(m,2) < zstart(2)-1 .or. bound_info%recids_c_loc(m,2) > zend(2)+1) then
bound_info%lcomprec_loc(m) = .true.
p = p + 1
end if
end if
bound_info%lskipsec_loc(m) = lbuf2(nn)
end do
end do
if (p > 0) then
write(*,*) "WARNING initibmwallfun: MPI rank", myid, "overrode", p, &
"facet section(s) to simple reconstruction because reconstruction cell falls outside halo range."
end if
deallocate(bound_info%bndpts)
if (myid == 0) then
deallocate(bound_info%secfacids)
deallocate(bound_info%secbndptids)
deallocate(bound_info%bnddst)
deallocate(bound_info%recpts)
deallocate(bound_info%recids_u)
deallocate(bound_info%recids_v)
deallocate(bound_info%recids_w)
deallocate(bound_info%recids_c)
deallocate(bound_info%lcomprec)
deallocate(bound_info%lskipsec)
end if
deallocate(lbndptsrank)
deallocate(ibuf1, ibuf2, rbuf1, rbuf2, rbuf3, ibufu, ibufv, ibufw, ibufc, lbuf1, lbuf2)
end subroutine initibmwallfun