subroutine read_forest(file_handle)
use mod_forest
use mod_global_parameters
use mod_space_filling_curve
use mod_amr_solution_node, only: getnode
use mod_functions_connectivity, only: get_level_range,build_connectivity,&
getigrids
use mod_amr_neighbors, only: asign_tree_neighbor
integer, intent(in) :: file_handle
integer, dimension(MPI_STATUS_SIZE) :: status
!integer :: ig^D, level, size_logical, Morton_no, igrid, ipe
integer :: ig1,ig2,ig3, level, Morton_no, igrid, ipe, isfc
logical :: leafbuf(nleafbuf)
integer :: ibuf, nvalid, nread
Morton_no=0
ipe=0
level=1
nleafs_level(1:nlevelshi) = 0
nparents = 0
! force a fill on the first node
ibuf = 0
nvalid = 0
nread = 0
call get_Morton_range
call level1_Morton_order
do isfc=1,nglev1
ig1=sfc_iglevel1(1,isfc)
ig2=sfc_iglevel1(2,isfc)
ig3=sfc_iglevel1(3,isfc)
allocate(tree_root(ig1,ig2,ig3)%node)
nullify(tree_root(ig1,ig2,ig3)%node%parent%node)
call read_node(tree_root(ig1,ig2,ig3),ig1,ig2,ig3,level)
end do
! Hand back the read-ahead: leave rank 0's file pointer exactly where the
! last consumed node ends, which is where it sat when this routine read one
! node at a time. read_snapshot_old relies on the sequential position.
if (mype==0 .and. nread>ibuf) then
call MPI_FILE_SEEK(file_handle, -int(nread-ibuf,kind=MPI_OFFSET_KIND)* &
int(size_logical,kind=MPI_OFFSET_KIND), MPI_SEEK_CUR, ierrmpi)
end if
call get_level_range
! Rebuild tree connectivity
call getigrids
call build_connectivity
! This call is here to ensure the sfc array is initialized
call amr_Morton_order()
contains
recursive subroutine read_node(tree,ig1,ig2,ig3,level)
implicit none
type(tree_node_ptr) :: tree
integer, intent(in) :: ig1,ig2,ig3, level
logical :: leaf
integer :: ic1,ic2,ic3, child_ig1,child_ig2,child_ig3, child_level
if (ibuf == nvalid) call refill_leafbuf
ibuf = ibuf + 1
leaf = leafbuf(ibuf)
tree%node%leaf=leaf
tree%node%ig1=ig1;tree%node%ig2=ig2;tree%node%ig3=ig3;
tree%node%level=level
tree%node%active=.true. .and. leaf
do ic3=1,2
do ic2=1,2
do ic1=1,2
nullify(tree%node%child(ic1,ic2,ic3)%node)
end do
end do
end do
nullify(tree%node%neighbor(1,1)%node,tree%node%neighbor(1,2)%node,&
tree%node%neighbor(1,3)%node,tree%node%neighbor(2,1)%node,&
tree%node%neighbor(2,2)%node,tree%node%neighbor(2,3)%node)
nullify(tree%node%next%node,tree%node%prev%node)
call asign_tree_neighbor(tree)
if (leaf) then
call add_to_linked_list(level,tree)
nleafs_level(level) = nleafs_level(level) + 1
Morton_no=Morton_no+1
if (Morton_no>Morton_stop(ipe)) ipe=ipe+1
igrid=getnode(ipe)
tree%node%igrid=igrid
tree%node%ipe=ipe
igrid_to_node(igrid,ipe)%node => tree%node
if (ipe==mype) sfc_to_igrid(Morton_no)=igrid
else
nparents = nparents + 1
tree%node%igrid=0
tree%node%ipe=-1
child_level=level+1
do ic3=1,2
do ic2=1,2
do ic1=1,2
child_ig1=2*(ig1-1)+ic1;child_ig2=2*(ig2-1)+ic2
child_ig3=2*(ig3-1)+ic3;
allocate(tree%node%child(ic1,ic2,ic3)%node)
tree%node%child(ic1,ic2,ic3)%node%parent%node => tree%node
call read_node(tree%node%child(ic1,ic2,ic3),child_ig1,child_ig2,&
child_ig3,child_level)
end do
end do
end do
end if
end subroutine read_node
!> Pull the next chunk of leaf flags in, one read and one broadcast for
!> up to nleafbuf nodes instead of one of each per node. Every rank walks
!> the tree in lockstep - the recursion is driven only by flags they all
!> hold - so they run out of buffer at the same node and this stays
!> collective without any agreement beyond that.
subroutine refill_leafbuf
implicit none
if (mype==0) then
call mpi_file_read_wrapper(file_handle,leafbuf,nleafbuf,MPI_LOGICAL,&
status,ierrmpi)
call MPI_GET_COUNT(status,MPI_LOGICAL,nread,ierrmpi)
end if
if (npe>1) call MPI_BCAST(leafbuf,nleafbuf,MPI_LOGICAL,0,icomm,ierrmpi)
! A short read at end of file is fine: the walk stops when the tree is
! complete and never looks past the last node the file holds, so the
! tail of the buffer is not consumed. nvalid therefore counts the whole
! buffer, and only rank 0 tracks what it actually read, for the rewind.
ibuf = 0
nvalid = nleafbuf
end subroutine refill_leafbuf
end subroutine read_forest