subroutine initibmwallfun(fname_bnd, fname_sec, dir, bound_info)
use modglobal, only : ifinput, ib, ie, itot, ih, jb, je, jtot, jh, kb, 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
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%nbndpts) :: lbndptsrank
logical, dimension(bound_info%nfctsecs) :: lfctsecsrank
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
allocate(bound_info%bndpts(bound_info%nbndpts,3))
! read u points
if (myid == 0) then
open (ifinput, file=fname_bnd)
read (ifinput, '(a80)') chmess
do n = 1, bound_info%nbndpts
read (ifinput, *) bound_info%bndpts(n,1), bound_info%bndpts(n,2), bound_info%bndpts(n,3)
end do
close (ifinput)
end if
call MPI_BCAST(bound_info%bndpts, bound_info%nbndpts*3, MPI_INTEGER, 0, comm3d, mpierr)
! Determine whether points are on this rank
bound_info%nbndptsrank = 0
do n = 1, bound_info%nbndpts
if ((bound_info%bndpts(n,1) >= zstart(1) .and. bound_info%bndpts(n,1) <= zend(1)) .and. &
(bound_info%bndpts(n,2) >= zstart(2) .and. bound_info%bndpts(n,2) <= zend(2))) then
lbndptsrank(n) = .true.
bound_info%nbndptsrank = bound_info%nbndptsrank + 1
else
lbndptsrank(n) = .false.
end if
end do
!write(*,*) "rank ", myid, " has ", bound_info%nbndptsrank, "points from ", fname_bnd
! Store indices of points on current rank - only loop through these points
allocate(bound_info%bndptsrank(bound_info%nbndptsrank)) ! index in global list
allocate(bound_info%bndpts_loc(bound_info%nbndptsrank,3)) ! location
m = 0
do n = 1, bound_info%nbndpts
if (lbndptsrank(n)) then
i = bound_info%bndpts(n,1) - zstart(1) + 1
j = bound_info%bndpts(n,2) - zstart(2) + 1
k = bound_info%bndpts(n,3) - zstart(3) + 1
if ((i < ib) .or. (i > ie) .or. (j < jb) .or. (j > je)) then
write(*,*) "problem in initibmwallfun", i, j
stop 1
end if
m = m + 1
bound_info%bndptsrank(m) = n
bound_info%bndpts_loc(m,:) = (/bound_info%bndpts(n,1),bound_info%bndpts(n,2),bound_info%bndpts(n,3)/)
end if
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(m,1) == 0 .or. bound_info%recids_u(m,2) == 0 .or. bound_info%recids_u(m,3) == 0) then
bound_info%lskipsec(n) = .true.
cycle
end if
if (bound_info%recids_v(m,1) == 0 .or. bound_info%recids_v(m,2) == 0 .or. bound_info%recids_v(m,3) == 0) then
bound_info%lskipsec(n) = .true.
cycle
end if
if (bound_info%recids_w(m,1) == 0 .or. bound_info%recids_w(m,2) == 0 .or. bound_info%recids_w(m,3) == 0) then
bound_info%lskipsec(n) = .true.
cycle
end if
if (bound_info%recids_c(m,1) == 0 .or. bound_info%recids_c(m,2) == 0 .or. bound_info%recids_c(m,3) == 0) then
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
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)
bound_info%lskipsec_loc(m) = bound_info%lskipsec(n)
end if
end do
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)
end subroutine initibmwallfun