diff --git a/src/boundaries.F90 b/src/boundaries.F90 index f66052e..7a28a95 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -3375,6 +3375,585 @@ module boundaries ! !=============================================================================== ! +! subroutine BOUNDARIES_EDGE_RESTRICT: +! ----------------------------------- +! +! Subroutine scans over all leaf blocks in order to find edge neighbors at +! the different levels, and perform the update of the edge boundaries between +! them. +! +! Arguments: +! +! idir - the direction to be processed; +! +!=============================================================================== +! + subroutine boundaries_edge_restrict(idir) + +! import external procedures and variables +! + use blocks , only : nsides + use blocks , only : block_meta, block_data + use blocks , only : list_meta + use blocks , only : block_info, pointer_info + use coordinates , only : ng + use coordinates , only : in , jn , kn + use coordinates , only : im , jm , km + use coordinates , only : ib , jb , kb + use coordinates , only : ie , je , ke + use coordinates , only : ibl, jbl, kbl + use coordinates , only : ieu, jeu, keu +#ifdef MPI + use equations , only : nv +#endif /* MPI */ + use mpitools , only : nproc, nprocs, npmax +#ifdef MPI + use mpitools , only : send_real_array, receive_real_array +#endif /* MPI */ + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer, intent(in) :: idir + +! local pointers +! + type(block_meta), pointer :: pmeta, pneigh +#ifdef MPI + type(block_info), pointer :: pinfo +#endif /* MPI */ + +! local variables +! + integer :: i , j , k + integer :: ih, jh, kh + integer :: il, jl, kl + integer :: iu, ju, ku + integer :: iret +#ifdef MPI + integer :: isend, irecv, nblocks, itag, l + +! local pointer arrays +! + type(pointer_info), dimension(0:npmax,0:npmax) :: block_array + +! local arrays +! + integer , dimension(0:npmax,0:npmax) :: block_counter + real(kind=8), dimension(:,:,:,:,:) , allocatable :: rbuf +#endif /* MPI */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for restrict boundary update +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate half sizes +! + ih = in / 2 + jh = jn / 2 +#if NDIMS == 3 + kh = kn / 2 +#endif /* NDIMS == 3 */ + +#ifdef MPI +!! 1. PREPARE THE BLOCK EXCHANGE ARRAYS FOR MPI +!! +! reset the exchange block counters +! + block_counter(:,:) = 0 + +! nullify the info pointers +! + do irecv = 0, npmax + do isend = 0, npmax + nullify(block_array(isend,irecv)%ptr) + end do + end do +#endif /* MPI */ + +!! 2. UPDATE VARIABLE EDGE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT +!! PROCESS AND PREPARE THE EXCHANGE BLOCK LIST OF BLOCKS WHICH BELONG TO +!! DIFFERENT PROCESSES +!! +! assign the pointer to the first block on the meta block list +! + pmeta => list_meta + +! scan all meta blocks and process blocks at the current level +! + do while(associated(pmeta)) + +! check if the block is leaf +! + if (pmeta%leaf) then + +! scan over all block corners +! +#if NDIMS == 3 + do k = 1, nsides +#endif /* NDIMS == 3 */ + do j = 1, nsides + do i = 1, nsides + +! assign pneigh to the current neighbor +! +#if NDIMS == 2 + pneigh => pmeta%edges(i,j,idir)%ptr +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + pneigh => pmeta%edges(i,j,k,idir)%ptr +#endif /* NDIMS == 3 */ + +! check if the neighbor is associated +! + if (associated(pneigh)) then + +! check if the neighbor is at higher level +! + if (pneigh%level > pmeta%level) then + +! skip if the block and its neighbor are not marked for update +! + if (pmeta%update .and. 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 (pmeta%process == nproc) then +#endif /* MPI */ + +! prepare the region indices for edge boundary update +! + if (i == 1) then + il = 1 + iu = ibl + else + il = ieu + iu = im + end if + if (j == 1) then + jl = 1 + ju = jbl + else + jl = jeu + ju = jm + end if +#if NDIMS == 3 + if (k == 1) then + kl = 1 + ku = kbl + else + kl = keu + ku = km + end if +#endif /* NDIMS == 3 */ + +! extract the corresponding edge region from the neighbor and insert it in +! the current data block +! + select case(idir) + case(1) + if (i == 1) then + il = ib + iu = ib + ih - 1 + else + il = ie - ih + 1 + iu = ie + end if + case(2) + if (j == 1) then + jl = jb + ju = jb + jh - 1 + else + jl = je - jh + 1 + ju = je + end if +#if NDIMS == 3 + case(3) + if (k == 1) then + kl = kb + ku = kb + kh - 1 + else + kl = ke - kh + 1 + ku = ke + end if +#endif /* NDIMS == 3 */ + end select +#if NDIMS == 2 + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv, 1:im, 1:jm, 1:km) & + , pmeta%data%q(1:nv,il:iu,jl:ju, 1:km)) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv, 1:im, 1:jm, 1:km) & + , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) +#endif /* NDIMS == 3 */ + +#ifdef MPI + end if ! pneigh on the current process + + else ! block and neighbor belong to different processes + +! increase the counter for number of blocks to exchange +! + block_counter(pneigh%process,pmeta%process) = & + block_counter(pneigh%process,pmeta%process) + 1 + +! allocate a new info object +! + allocate(pinfo) + +! fill out only fields which are used +! + pinfo%block => pmeta + pinfo%neigh => pneigh + pinfo%direction = idir + pinfo%corner(1) = i + pinfo%corner(2) = j +#if NDIMS == 3 + pinfo%corner(3) = k +#endif /* NDIMS == 3 */ + +! nullify pointer fields of the object +! + nullify(pinfo%prev) + nullify(pinfo%next) + +! if the list is not empty append the newly created block to it +! + if (associated(block_array(pneigh%process & + ,pmeta%process)%ptr)) & + pinfo%prev => block_array(pneigh%process & + ,pmeta%process)%ptr + +! point the list to the newly created block +! + block_array(pneigh%process,pmeta%process)%ptr => pinfo + + end if ! block and neighbor belong to different processes +#endif /* MPI */ + + end if ! pmeta and pneigh marked for update + + end if ! neighbor at the same level + + end if ! neighbor associated + + end do ! i = 1, nsides + end do ! j = 1, nsides +#if NDIMS == 3 + end do ! k = 1, nsides +#endif /* NDIMS == 3 */ + + end if ! leaf + +! associate the pointer to the next meta block +! + pmeta => pmeta%next + + end do ! meta blocks + +#ifdef MPI +!! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES +!! +! iterate over sending and receiving processors +! + do irecv = 0, npmax + do isend = 0, npmax + +! process only pairs which have something to exchange +! + if (block_counter(isend,irecv) > 0) then + +! obtain the number of blocks to exchange +! + nblocks = block_counter(isend,irecv) + +! prepare the tag for communication +! + itag = 100 * (irecv * nprocs + isend + 1) + 22 + +! allocate data buffer for variables to exchange +! + select case(idir) +#if NDIMS == 2 + case(1) + allocate(rbuf(nblocks,nv,ih,ng,km)) + case(2) + allocate(rbuf(nblocks,nv,ng,jh,km)) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + case(1) + allocate(rbuf(nblocks,nv,ih,ng,ng)) + case(2) + allocate(rbuf(nblocks,nv,ng,jh,ng)) + case(3) + allocate(rbuf(nblocks,nv,ng,ng,kh)) +#endif /* NDIMS == 3 */ + end select + +! if isend == nproc we are sending data from the neighbor block +! + if (isend == nproc) then + +! reset the block counter +! + l = 0 + +! associate the pointer with the first block in the exchange list +! + pinfo => block_array(isend,irecv)%ptr + +! scan over all blocks on the block exchange list +! + do while(associated(pinfo)) + +! increase the block counter +! + l = l + 1 + +! assign pneigh to the associated neighbor block +! + pneigh => pinfo%neigh + +! get the corner coordinates +! + i = pinfo%corner(1) + j = pinfo%corner(2) +#if NDIMS == 3 + k = pinfo%corner(3) +#endif /* NDIMS == 3 */ + +! extract the corresponding edge region from the neighbor and insert it +! to the buffer +! + select case(idir) + case(1) +#if NDIMS == 2 + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv,1:im,1:jm,1:km) & + , rbuf(l,1:nv,1:ih,1:ng,1:km)) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv,1:im,1:jm,1:km) & + , rbuf(l,1:nv,1:ih,1:ng,1:ng)) +#endif /* NDIMS == 3 */ + case(2) +#if NDIMS == 2 + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv,1:im,1:jm,1:km) & + , rbuf(l,1:nv,1:ng,1:jh,1:km)) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv,1:im,1:jm,1:km) & + , rbuf(l,1:nv,1:ng,1:jh,1:ng)) +#endif /* NDIMS == 3 */ +#if NDIMS == 3 + case(3) + call block_edge_restrict(idir, i, j, k & + , pneigh%data%q(1:nv,1:im,1:jm,1:km) & + , rbuf(l,1:nv,1:ng,1:ng,1:kh)) +#endif /* NDIMS == 3 */ + end select + +! associate the pointer with the next block +! + pinfo => pinfo%prev + + end do ! %ptr block list + +! send the data buffer to another process +! + call send_real_array(size(rbuf), irecv, itag, rbuf(:,:,:,:,:), iret) + + end if ! isend = nproc + +! if irecv == nproc we are receiving data from the neighbor block +! + if (irecv == nproc) then + +! receive the data buffer +! + call receive_real_array(size(rbuf(:,:,:,:,:)), isend, itag & + , rbuf(:,:,:,:,:), iret) + +! reset the block counter +! + l = 0 + +! associate the pointer with the first block in the exchange list +! + pinfo => block_array(isend,irecv)%ptr + +! iterate over all received blocks and update boundaries of the corresponding +! data blocks +! + do while(associated(pinfo)) + +! increase the block counter +! + l = l + 1 + +! assign a pointer to the associated data block +! + pmeta => pinfo%block + +! get the corner coordinates +! + i = pinfo%corner(1) + j = pinfo%corner(2) +#if NDIMS == 3 + k = pinfo%corner(3) +#endif /* NDIMS == 3 */ + +! calculate the insertion indices +! + if (i == 1) then + il = 1 + iu = ibl + else + il = ieu + iu = im + end if + if (j == 1) then + jl = 1 + ju = jbl + else + jl = jeu + ju = jm + end if +#if NDIMS == 3 + if (k == 1) then + kl = 1 + ku = kbl + else + kl = keu + ku = km + end if +#endif /* NDIMS == 3 */ + +! update the corresponding corner region of the current block +! + select case(idir) + case(1) + if (i == 1) then + il = ib + iu = ib + ih - 1 + else + il = ie - ih + 1 + iu = ie + end if +#if NDIMS == 2 + pmeta%data%q(1:nv,il:iu,jl:ju, 1:km) = & + rbuf(l,1:nv,1:ih,1:ng,1:km) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & + rbuf(l,1:nv,1:ih,1:ng,1:ng) +#endif /* NDIMS == 3 */ + case(2) + if (j == 1) then + jl = jb + ju = jb + jh - 1 + else + jl = je - jh + 1 + ju = je + end if +#if NDIMS == 2 + pmeta%data%q(1:nv,il:iu,jl:ju, 1:km) = & + rbuf(l,1:nv,1:ng,1:jh,1:km) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & + rbuf(l,1:nv,1:ng,1:jh,1:ng) +#endif /* NDIMS == 3 */ +#if NDIMS == 3 + case(3) + if (k == 1) then + kl = kb + ku = kb + kh - 1 + else + kl = ke - kh + 1 + ku = ke + end if + pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & + rbuf(l,1:nv,1:ng,1:ng,1:kh) +#endif /* NDIMS == 3 */ + end select + +! associate the pointer with the next block +! + pinfo => pinfo%prev + + end do ! %ptr block list + + end if ! irecv = nproc + +! deallocate data buffer +! + if (allocated(rbuf)) deallocate(rbuf) + +! associate the pointer with the first block in the exchange list +! + pinfo => block_array(isend,irecv)%ptr + +! scan over all blocks on the exchange block list +! + do while(associated(pinfo)) + +! associate the exchange list pointer +! + block_array(isend,irecv)%ptr => pinfo%prev + +! nullify the pointer fields +! + nullify(pinfo%prev) + nullify(pinfo%next) + nullify(pinfo%block) + nullify(pinfo%neigh) + +! deallocate the object +! + deallocate(pinfo) + +! associate the pointer with the next block +! + pinfo => block_array(isend,irecv)%ptr + + end do ! %ptr block list + + end if ! if block_count > 0 + + end do ! isend + end do ! irecv +#endif /* MPI */ + +#ifdef PROFILE +! stop accounting time for restrict boundary update +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine boundaries_edge_restrict +! +!=============================================================================== +! ! subroutine BOUNDARIES_CORNER_COPY: ! --------------------------------- !