subroutine initibmwallfun(fname_bnd, fname_sec, dir, bound_info)
use modglobal, only : ifinput, ib, ie, itot, ih, jb, je, jtot, jh, kb, ke, ktot, kh, &
xf, yf, zf, xh, yh, zh, dx, dy, dzh, 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
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, checkxl, checkxu, checkyl, checkyu, checkzl, checkzu, inter_dist
integer i, j, k, n, m, norm_align, dir_align, pos, p
real dst
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.
do m = 1, bound_info%nbndptsrank
lbndptsrank(bound_info%bndptsrank(m)) = .true.
end do
allocate(bound_info%secfacids(bound_info%nfctsecs))
allocate(bound_info%secareas(bound_info%nfctsecs))
allocate(bound_info%secbndptids(bound_info%nfctsecs))
!allocate(bound_info%intpts(bound_info%nfctsecs,3))
allocate(bound_info%bnddst(bound_info%nfctsecs))
!allocate(bound_info%bndvec(bound_info%nfctsecs,3))
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))
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
call MPI_BCAST(bound_info%secfacids, bound_info%nfctsecs, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%secareas, bound_info%nfctsecs, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%secbndptids, bound_info%nfctsecs, MPI_INTEGER, 0, comm3d, mpierr)
!call MPI_BCAST(bound_info%intpts, bound_info%nfctsecs*3, MY_REAL, 0, comm3d, mpierr)
!call MPI_BCAST(bound_info%bndvec, bound_info%nfctsecs*3, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%bnddst, bound_info%nfctsecs, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%recpts, bound_info%nfctsecs*3, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%recids_u, bound_info%nfctsecs*3, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%recids_v, bound_info%nfctsecs*3, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%recids_w, bound_info%nfctsecs*3, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%recids_c, bound_info%nfctsecs*3, MPI_INTEGER, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%lskipsec, bound_info%nfctsecs, MPI_LOGICAL, 0, comm3d, mpierr)
call MPI_BCAST(bound_info%lcomprec, bound_info%nfctsecs, MPI_LOGICAL, 0, comm3d, mpierr)
! Determine whether section needs to be updated by this rank
bound_info%nfctsecsrank = 0
do n = 1, bound_info%nfctsecs
if (lbndptsrank(bound_info%secbndptids(n))) then
lfctsecsrank(n) = .true.
bound_info%nfctsecsrank = bound_info%nfctsecsrank + 1
else
lfctsecsrank(n) = .false.
end if
end do
! Store indices of sections on current rank - only loop through these sections
allocate(bound_info%fctsecsrank(bound_info%nfctsecsrank))
! allocate local arrays
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))
m = 0
p = 0 ! counter for sections overridden to simple reconstruction
do n = 1, bound_info%nfctsecs
if (lfctsecsrank(n)) then
m = m + 1
bound_info%fctsecsrank(m) = n
bound_info%secfacids_loc(m) = bound_info%secfacids(n) ! facet id
bound_info%secareas_loc(m) = bound_info%secareas(n)
bound_info%secbndpts_loc(m,:) = bound_info%bndpts(bound_info%secbndptids(n),:) ! boundary point location (in global coordinates)
if (bound_info%bndpts(bound_info%secbndptids(n),1) < zstart(1) .or. bound_info%bndpts(bound_info%secbndptids(n),1) > zend(1)) then
write(*,*) "problem in x boundary points on : ", myid, n, bound_info%secbndptids(n), bound_info%bndpts(bound_info%secbndptids(n),1), zstart(1), zend(1)
end if
if (bound_info%bndpts(bound_info%secbndptids(n),2) < zstart(2) .or. bound_info%bndpts(bound_info%secbndptids(n),2) > zend(2)) then
write(*,*) "problem in y boundary points on rank: ", myid, n, bound_info%secbndptids(n), bound_info%bndpts(bound_info%secbndptids(n),2), zstart(2), zend(2)
end if
bound_info%bnddst_loc(m) = bound_info%bnddst(n)
bound_info%recpts_loc(m,:) = bound_info%recpts(n,:)
bound_info%recids_u_loc(m,:) = bound_info%recids_u(n,:)
bound_info%recids_v_loc(m,:) = bound_info%recids_v(n,:)
bound_info%recids_w_loc(m,:) = bound_info%recids_w(n,:)
bound_info%recids_c_loc(m,:) = bound_info%recids_c(n,:)
bound_info%lcomprec_loc(m) = bound_info%lcomprec(n)
! If any reconstruction cell index falls outside this rank's halo-accessible range,
! override to simple reconstruction to avoid out-of-bounds access in trilinear_interp_var.
!
! Skipped sections (lskipsec(n)) never reach trilinear_interp_var
! so they are harmless even if their reconstruction cell is out of range.
!
! These sections then self-skip at runtime via the log(dist/z0) <= 1 check, contributing zero stress — a small, localized inaccuracy
if (.not. bound_info%lcomprec_loc(m) .and. .not. bound_info%lskipsec(n)) 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) = bound_info%lskipsec(n)
end if
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)
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)
deallocate(lbndptsrank)
end subroutine initibmwallfun