BOUNDARIES: Rewrite boundaries_face_restrict().

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-02-05 18:05:11 -03:00
parent 65be245c97
commit cdcf90ae30

View File

@ -912,7 +912,7 @@ module boundaries
! restrict face boundaries from higher level blocks
!
do idir = 1, ndims
call boundaries_face_restrict(idir)
call boundaries_face_restrict(idir, status)
end do ! idir
#endif /* NDIMS == 3 */
@ -1489,180 +1489,175 @@ module boundaries
! subroutine BOUNDARIES_FACE_RESTRICT:
! -----------------------------------
!
! Subroutine updates the face boundaries from blocks on higher level.
! Subroutine updates the blocks' face ghost zones from blocks at
! higher levels.
!
! Arguments:
!
! idir - the direction to be processed;
! idir - the direction of the update;
! status - the call status;
!
!===============================================================================
!
subroutine boundaries_face_restrict(idir)
subroutine boundaries_face_restrict(idir, status)
! import external procedures and variables
!
use blocks , only : nsides
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
use coordinates, only : bcells, ncells_half, nb, ne
#ifdef MPI
use coordinates, only : nh => ncells_half, ng => nghosts
use coordinates, only : nghosts
#endif /* MPI */
use coordinates, only : faces_gr
use equations , only : nv
use helpers , only : print_message
#ifdef MPI
use mpitools , only : nproc, npairs, pairs
use mpitools , only : exchange_arrays
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(in) :: idir
integer, intent(in) :: idir
integer, intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta, pneigh
type(block_data), pointer :: pdata
type(block_leaf), pointer :: pleaf
#ifdef MPI
type(block_info), pointer :: pinfo
#endif /* MPI */
! local variables
!
integer :: i , j , k
integer :: il, jl, kl
integer :: iu, ju, ku
integer :: i, il, iu
integer :: j, jl, ju
integer :: k, kl, ku
#ifdef MPI
integer :: sproc, rproc
integer :: scount, rcount
integer :: l, p
! local arrays
!
integer, dimension(3) :: dm, pm
real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf
#endif /* MPI */
!
logical , save :: first = .true.
integer, dimension(2,2), save :: nlims, tlims
character(len=*), parameter :: loc = &
'BOUNDARIES::boundaries_face_restrict()'
!-------------------------------------------------------------------------------
!
if (first) then
nlims(1,1) = 1
nlims(1,2) = nb - 1
nlims(2,1) = ne + 1
nlims(2,2) = bcells
tlims(1,1) = nb
tlims(1,2) = nb + ncells_half - 1
tlims(2,1) = ne - ncells_half + 1
tlims(2,2) = ne
first = .false.
end if
status = 0
#ifdef MPI
sproc = 0
rproc = 0
scount = 0
rcount = 0
! prepare the block exchange structures
!
dm( : ) = ncells_half
dm(idir) = nghosts
call prepare_exchange_array()
#endif /* MPI */
! update boundaries between blocks on the same process
!
! associate pleaf with the first block on the leaf list
!
pleaf => list_leaf
! scan all leaf meta blocks in the list
!
do while(associated(pleaf))
! get the associated meta block
!
pmeta => pleaf%meta
pdata => pmeta%data
! scan over all block corners
!
do k = 1, nsides
pm(3) = k
if (idir == 3) then
kl = nlims(pm(3),1)
ku = nlims(pm(3),2)
else
kl = tlims(pm(3),1)
ku = tlims(pm(3),2)
end if
do j = 1, nsides
pm(2) = j
if (idir == 2) then
jl = nlims(pm(2),1)
ju = nlims(pm(2),2)
else
jl = tlims(pm(2),1)
ju = tlims(pm(2),2)
end if
do i = 1, nsides
pm(1) = i
if (idir == 1) then
il = nlims(pm(1),1)
iu = nlims(pm(1),2)
else
il = tlims(pm(1),1)
iu = tlims(pm(1),2)
end if
! associate pneigh with the current neighbor
!
pneigh => pmeta%faces(i,j,k,idir)%ptr
! check if the neighbor is associated
!
if (associated(pneigh)) then
! check if the neighbor is at higher level
!
if (pneigh%level > pmeta%level) then
! process only blocks and neighbors which are marked for update
!
if (pmeta%update .or. pneigh%update) then
#ifdef MPI
! check if the block and its neighbor belong to the same process
!
if (pmeta%process == pneigh%process) then
! check if the neighbor belongs to the current process
!
if (pneigh%process == nproc) then
#endif /* MPI */
! prepare region indices of the block and its neighbor for the face boundary
! update
!
il = faces_gr(i,j,k,idir)%l(1)
jl = faces_gr(i,j,k,idir)%l(2)
kl = faces_gr(i,j,k,idir)%l(3)
iu = faces_gr(i,j,k,idir)%u(1)
ju = faces_gr(i,j,k,idir)%u(2)
ku = faces_gr(i,j,k,idir)%u(3)
! extract the corresponding face region from the neighbor and insert it in
! the current data block
!
call block_face_restrict(idir, [ i, j, k ] &
, pneigh%data%q(1:nv, : , : , : ) &
, pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku))
call block_face_restrict(idir, pm, pneigh%data%q, &
pdata%q(:,il:iu,jl:ju,kl:ku))
#ifdef MPI
end if ! pneigh on the current process
end if
else ! block and neighbor belong to different processes
else
! append the block to the exchange list
!
call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /))
call append_exchange_block(pmeta, pneigh, idir, pm)
end if ! block and neighbor belong to different processes
end if
#endif /* MPI */
end if ! pmeta and pneigh marked for update
end if
end if ! neighbor at the same level
end if
end if ! neighbor associated
end if
end do ! i = 1, nsides
end do ! j = 1, nsides
end do ! k = 1, nsides
end do
end do
end do
! associate pleaf with the next leaf on the list
!
pleaf => pleaf%next
end do ! over leaf blocks
end do
#ifdef MPI
!! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES
!!
! iterate over all process pairs
!
do p = 1, npairs
! process only pairs related to this process
!
if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then
! get sending and receiving process identifiers (depending on pair member)
!
if (pairs(p,1) == nproc) then
sproc = pairs(p,1)
rproc = pairs(p,2)
@ -1672,157 +1667,88 @@ module boundaries
rproc = pairs(p,1)
end if
! get the number of blocks to exchange
!
scount = bcount(sproc,rproc)
rcount = bcount(rproc,sproc)
! process only pairs which have anything to exchange
!
if ((scount + rcount) > 0) then
if (scount > 0 .or. rcount > 0) then
! allocate data buffer for variables to exchange
!
select case(idir)
case(1)
allocate(sbuf(scount,nv,ng,nh,nh))
allocate(rbuf(rcount,nv,ng,nh,nh))
case(2)
allocate(sbuf(scount,nv,nh,ng,nh))
allocate(rbuf(rcount,nv,nh,ng,nh))
case(3)
allocate(sbuf(scount,nv,nh,nh,ng))
allocate(rbuf(rcount,nv,nh,nh,ng))
end select
allocate(sbuf(nv,dm(1),dm(2),dm(3),scount), &
rbuf(nv,dm(1),dm(2),dm(3),rcount), stat=status)
if (status /= 0) &
call print_message(loc, "Could not allocate the exchange buffers!")
!! PREPARE BLOCKS FOR SENDING
!!
! reset the block counter
!
l = 0
if (scount > 0) then
l = 0
! associate the pointer with the first block in the exchange list
!
pinfo => barray(sproc,rproc)%ptr
pinfo => barray(sproc,rproc)%ptr
! scan over all blocks on the block exchange list
!
do while(associated(pinfo))
do while(associated(pinfo))
! increase the block counter
!
l = l + 1
l = l + 1
! assign pneigh to the associated neighbor block
!
pneigh => pinfo%neigh
pdata => pinfo%neigh%data
! get the corner coordinates
!
i = pinfo%corner(1)
j = pinfo%corner(2)
k = pinfo%corner(3)
call block_face_restrict(idir, pinfo%corner, &
pdata%q, sbuf(:,:,:,:,l))
! restrict the corresponding face region from the neighbor and insert it
! to the buffer
!
select case(idir)
case(1)
call block_face_restrict(idir, [ i, j, k ] &
, pneigh%data%q(1:nv, : , : , : ) &
, sbuf(l,1:nv,1:ng,1:nh,1:nh))
case(2)
call block_face_restrict(idir, [ i, j, k ] &
, pneigh%data%q(1:nv, : , : , : ) &
, sbuf(l,1:nv,1:nh,1:ng,1:nh))
case(3)
call block_face_restrict(idir, [ i, j, k ] &
, pneigh%data%q(1:nv, : , : , : ) &
, sbuf(l,1:nv,1:nh,1:nh,1:ng))
end select
pinfo => pinfo%prev
end do
end if
! associate the pointer with the next block
!
pinfo => pinfo%prev
end do ! %ptr block list
!! SEND PREPARED BLOCKS AND RECEIVE NEW ONES
!!
! exchange data
!
call exchange_arrays(rproc, p, sbuf, rbuf)
!! PROCESS RECEIVED BLOCKS
!!
! reset the block counter
!
l = 0
if (rcount > 0) then
l = 0
! associate the pointer with the first block in the exchange list
!
pinfo => barray(rproc,sproc)%ptr
pinfo => barray(rproc,sproc)%ptr
! iterate over all received blocks and update boundaries of the corresponding
! data blocks
!
do while(associated(pinfo))
do while(associated(pinfo))
! increase the block counter
!
l = l + 1
l = l + 1
! assign a pointer to the associated data block
!
pmeta => pinfo%meta
pmeta => pinfo%meta
! get the corner coordinates
!
i = pinfo%corner(1)
j = pinfo%corner(2)
k = pinfo%corner(3)
pm(:) = pinfo%corner(:)
! prepare region indices for the face boundary update
!
il = faces_gr(i,j,k,idir)%l(1)
jl = faces_gr(i,j,k,idir)%l(2)
kl = faces_gr(i,j,k,idir)%l(3)
iu = faces_gr(i,j,k,idir)%u(1)
ju = faces_gr(i,j,k,idir)%u(2)
ku = faces_gr(i,j,k,idir)%u(3)
if (idir == 1) then
il = nlims(pm(1),1)
iu = nlims(pm(1),2)
else
il = tlims(pm(1),1)
iu = tlims(pm(1),2)
end if
if (idir == 2) then
jl = nlims(pm(2),1)
ju = nlims(pm(2),2)
else
jl = tlims(pm(2),1)
ju = tlims(pm(2),2)
end if
if (idir == 3) then
kl = nlims(pm(3),1)
ku = nlims(pm(3),2)
else
kl = tlims(pm(3),1)
ku = tlims(pm(3),2)
end if
! update the corresponding face region of the current block
!
select case(idir)
case(1)
pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = &
rbuf(l,1:nv,1:ng,1:nh,1:nh)
case(2)
pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = &
rbuf(l,1:nv,1:nh,1:ng,1:nh)
case(3)
pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = &
rbuf(l,1:nv,1:nh,1:nh,1:ng)
end select
pmeta%data%q(:,il:iu,jl:ju,kl:ku) = rbuf(:,:,:,:,l)
! associate the pointer with the next block
!
pinfo => pinfo%prev
pinfo => pinfo%prev
end do
end if
end do ! %ptr block list
deallocate(sbuf, rbuf, stat=status)
if (status /= 0) &
call print_message(loc, &
"Could not deallocate the exchange buffers!")
! deallocate data buffer
!
deallocate(sbuf, rbuf)
end if
end if ! (scount + rcount) > 0
end if
end if ! pairs(p,1) == nproc || pairs(p,2) == nproc
end do
end do ! p = 1, npairs
! release the memory used by the array of exchange block lists
!
call release_exchange_array()
#endif /* MPI */
@ -2026,9 +1952,7 @@ module boundaries
pdata => pinfo%neigh%data
pm(:) = pinfo%corner(:)
call block_face_prolong(idir, pm, pdata%q, &
call block_face_prolong(idir, pinfo%corner, pdata%q, &
sbuf(:,:,:,:,l), status)
pinfo => pinfo%prev