forest.t

Go to the documentation of this file.

!> build root forest
subroutine init_forest_root
  use mod_forest
  use mod_global_parameters
  use mod_space_filling_curve
 
  integer :: ig^D, level, igrid, ipe
  integer :: iside, i^D, Morton_no, isfc
 
  integer, external :: getnode
 
  level=1
  morton_no=0
  ipe=0
  nparents=0
  nleafs={ng^d(1)*}
  nleafs_active=nleafs
  nleafs_level(1)={ng^d(1)*}
  nleafs_level(2:nlevelshi)=0
  call get_morton_range
  ! Generate Morton-order space-filling curve to connect level 1 blocks
  call level1_morton_order
  do isfc=1,nglev1
     ig^d=sfc_iglevel1(^d,isfc)\
     morton_no=morton_no+1
     if (morton_no>morton_stop(ipe)) ipe=ipe+1
     igrid=getnode(ipe)
     if (ipe==mype) then
        sfc_to_igrid(morton_no)=igrid
        igrid_to_sfc(igrid)=morton_no
     end if
     call init_tree_leaf(tree_root(ig^d),ig^d,level,igrid,ipe,.true.)
  end do
 
  ! update root neighbor
  {do ig^db=1,ng^db(1)\}
     {do iside=1,2
        i^dd=kr(^dd,^d)*(2*iside-3);
        call find_root_neighbor(tree_root(ig^dd)%node%neighbor(iside,^d), &
                                tree_root(ig^dd),i^dd)
     end do\}
  {end do\}
 
  ! This call is here to ensure the sfc array is initialized
  call amr_morton_order()
 
end subroutine init_forest_root
 
subroutine init_tree_leaf(tree,ig^D,level,igrid,ipe,active)
  use mod_forest
  implicit none
 
  type(tree_node_ptr) :: tree
  integer, intent(in) :: ig^D, level, igrid, ipe
  logical, intent(in) :: active
  integer :: ic^D
 
  allocate(tree%node)
 
  tree%node%ig^d=ig^d;
  tree%node%level=level
  tree%node%igrid=igrid
  tree%node%ipe=ipe
 
  tree%node%leaf=.true.
  tree%node%active=active
 
  nullify(tree%node%parent%node)
  {do ic^db=1,2\}
     nullify(tree%node%child(ic^d)%node)
  {end do\}
 
  call add_to_linked_list(level,tree)
 
  ! initialize neighbor pointers
  nullify({tree%node%neighbor(1,^d)%node},{tree%node%neighbor(2,^d)%node})
 
  igrid_to_node(igrid,ipe)%node => tree%node
 
end subroutine init_tree_leaf
 
subroutine coarsen_tree_leaf(igrid,ipe,child_igrid,child_ipe,active)
  use mod_forest
  implicit none
 
  integer, intent(in) :: igrid, ipe
  integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
  logical, intent(out) :: active
 
 
  integer :: level, ic^D, child_level, iside, iotherside, vote
  type(tree_node_ptr) :: tree, child, child_neighbor
 
  tree%node => igrid_to_node(child_igrid(1^d&),child_ipe(1^d&))%node%parent%node
  level=tree%node%level
 
  call add_to_linked_list(level,tree)
 
  child_level=level+1
  vote=0
 
  {do ic^db=1,2\}
     child%node => tree%node%child(ic^d)%node
 
  !  vote for active:
     if(child%node%active) vote=vote+1
 
     call delete_from_linked_list(child_level,child)
 
     ! update neighbor pointers
     {iside=ic^d
     child_neighbor%node => child%node%neighbor(iside,^d)%node
     if (associated(child_neighbor%node)) then
        if (child%node%ig^d==child_neighbor%node%ig^d) then ! pole
           nullify(child_neighbor%node%neighbor(iside,^d)%node)
        else
           iotherside=3-iside
           nullify(child_neighbor%node%neighbor(iotherside,^d)%node)
        end if
     end if\}
 
     nullify(tree%node%child(ic^d)%node)
     deallocate(igrid_to_node(child_igrid(ic^d),child_ipe(ic^d))%node)
  {end do\}
 
  tree%node%leaf=.true.
  tree%node%igrid=igrid
  tree%node%ipe=ipe
  igrid_to_node(igrid,ipe)%node => tree%node
 
  !  Count the vote and set active/passive state:
 
  if (vote /= 2**^nd) then 
  !if (vote == 0) then 
     tree%node%active = .false.
     nleafs_active = nleafs_active - vote
  else
     tree%node%active = .true.
     nleafs_active = nleafs_active - vote + 1
  end if
  active = tree%node%active
 
  nleafs=nleafs-2**^nd+1
  nparents=nparents-1
  nleafs_level(child_level)=nleafs_level(child_level)-2**^nd
  nleafs_level(level)=nleafs_level(level)+1
 
end subroutine coarsen_tree_leaf
 
subroutine refine_tree_leaf(child_igrid,child_ipe,igrid,ipe,active)
  use mod_forest
  use mod_global_parameters
 
  integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
  integer, intent(in) :: igrid, ipe
  logical, intent(out):: active
 
  integer :: ig^D, level, i^D, ic^D, child_ig^D, child_level, iside
  integer :: my_neighbor_type
  logical, dimension(ndim) :: pole
  type(tree_node_ptr) :: tree, child, my_neighbor
 
  tree%node => igrid_to_node(igrid,ipe)%node
  ig^d=tree%node%ig^d;
  level=tree%node%level
  active=tree%node%active
 
  tree%node%ipe=-1
  tree%node%igrid=0
  tree%node%leaf=.false.
  tree%node%active=.true.
 
  call delete_from_linked_list(level,tree)
 
  child_level=level+1
 
  {do ic^db=1,2\}
     child_ig^d=2*(ig^d-1)+ic^d;
     call init_tree_leaf(child,child_ig^d,child_level, &
                               child_igrid(ic^d),child_ipe(ic^d),active)
 
     igrid_to_node(child_igrid(ic^d),child_ipe(ic^d))%node => child%node
 
     tree%node%child(ic^d)%node => child%node
     child%node%parent%node => tree%node
  {end do\}
 
  ! update neighbor pointers
  {do ic^db=1,2\}
     child%node => tree%node%child(ic^d)%node
     {iside=ic^d
     i^dd=kr(^dd,^d)*(2*iside-3);
     call find_neighbor(my_neighbor,my_neighbor_type,child,i^dd,pole)
     select case (my_neighbor_type)
     case (neighbor_sibling, neighbor_fine)
        child%node%neighbor(iside,^d)%node => my_neighbor%node
        if (pole(^d)) then
           my_neighbor%node%neighbor(iside,^d)%node => child%node
        else
           my_neighbor%node%neighbor(3-iside,^d)%node => child%node
        end if
     case default
        nullify(child%node%neighbor(iside,^d)%node)
     end select
     child%node%neighbor(3-ic^d,^d)%node=>tree%node%child(3-ic^d^d%ic^dd)%node\}
  {end do\}
 
  nleafs=nleafs+2**^nd-1
  nparents=nparents+1
  nleafs_level(child_level)=nleafs_level(child_level)+2**^nd
  nleafs_level(level)=nleafs_level(level)-1
 
  if (active) nleafs_active = nleafs_active + 2**^nd-1
 
end subroutine refine_tree_leaf
 
subroutine change_ipe_tree_leaf(recv_igrid,recv_ipe,send_igrid,send_ipe)
  use mod_forest
  implicit none
 
  integer, intent(in) :: recv_igrid, recv_ipe, send_igrid, send_ipe
 
  type(tree_node_ptr) :: tree
 
  tree%node => igrid_to_node(send_igrid,send_ipe)%node
 
  tree%node%igrid=recv_igrid
  tree%node%ipe=recv_ipe
 
  nullify(igrid_to_node(send_igrid,send_ipe)%node)
  igrid_to_node(recv_igrid,recv_ipe)%node => tree%node
 
end subroutine change_ipe_tree_leaf
 
subroutine add_to_linked_list(level,tree)
  use mod_forest
  implicit none
 
  integer, intent(in) :: level
  type(tree_node_ptr) :: tree
 
  nullify(tree%node%next%node)
  if (associated(level_head(level)%node)) then
     tree%node%prev%node => level_tail(level)%node
     level_tail(level)%node%next%node => tree%node
     level_tail(level)%node => tree%node
  else
     level_head(level)%node => tree%node
     level_tail(level)%node => tree%node
     nullify(tree%node%prev%node)
  end if
 
end subroutine add_to_linked_list
 
subroutine delete_from_linked_list(level,tree)
  use mod_forest
  implicit none
 
  integer, intent(in) :: level
  type(tree_node_ptr) :: tree
 
  type(tree_node_ptr) :: next, prev
 
  prev%node => tree%node%prev%node
  next%node => tree%node%next%node
  if (associated(next%node).and.associated(prev%node)) then
     prev%node%next%node => next%node
     next%node%prev%node => prev%node
  else if (associated(prev%node)) then
     level_tail(level)%node => prev%node
     nullify(prev%node%next%node)
  else if (associated(next%node)) then
     level_head(level)%node => next%node
     nullify(next%node%prev%node)
  else
     nullify(level_head(level)%node)
     nullify(level_tail(level)%node)
  end if
 
end subroutine delete_from_linked_list
 
subroutine write_forest(file_handle)
  use mod_forest
  use mod_global_parameters
 
  integer, intent(in) :: file_handle
 
  integer, dimension(MPI_STATUS_SIZE) :: status
  integer :: ig^D,isfc
 
  do isfc=1,nglev1
     ig^d=sfc_iglevel1(^d,isfc)\
     call write_node(tree_root(ig^d))
  end do
 
  contains
 
    recursive subroutine write_node(tree)
      implicit none
 
      type(tree_node_ptr) :: tree
 
      integer :: ic^d
 
      call mpi_file_write(file_handle,tree%node%leaf,1,mpi_logical,status,ierrmpi)
 
      if (.not.tree%node%leaf) then
         {do ic^db=1,2\}
            call write_node(tree%node%child(ic^d))
         {end do\}
      end if
 
    end subroutine write_node
 
end subroutine write_forest
 
subroutine read_forest(file_handle)
  use mod_forest
  use mod_global_parameters
  use mod_space_filling_curve
 
  integer, intent(in) :: file_handle
 
  integer, dimension(MPI_STATUS_SIZE) :: status
  !integer :: ig^D, level, size_logical, Morton_no, igrid, ipe
  integer :: ig^D, level, Morton_no, igrid, ipe, isfc
  integer, external :: getnode
 
  morton_no=0
  ipe=0
  level=1
  nleafs_level(1:nlevelshi) = 0
  nparents = 0
 
  call get_morton_range
  call level1_morton_order
  do isfc=1,nglev1
     ig^d=sfc_iglevel1(^d,isfc)\
     allocate(tree_root(ig^d)%node)
     nullify(tree_root(ig^d)%node%parent%node)
     call read_node(tree_root(ig^d),ig^d,level)
  end do
 
  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,ig^D,level)
      implicit none
 
      type(tree_node_ptr) :: tree
      integer, intent(in) :: ig^d, level
 
      logical :: leaf
      integer :: ic^d, child_ig^d, child_level
 
      if (mype==0) then
        call mpi_file_read(file_handle,leaf,1,mpi_logical, &
                                  status,ierrmpi)
      end if
      if (npe>1)  call mpi_bcast(leaf,1,mpi_logical,0,icomm,ierrmpi)
 
      tree%node%leaf=leaf
      tree%node%ig^d=ig^d;
      tree%node%level=level
      tree%node%active=.true. .and. leaf
 
      {do ic^db=1,2\}
         nullify(tree%node%child(ic^d)%node)
      {end do\}
      nullify({tree%node%neighbor(1,^d)%node},{tree%node%neighbor(2,^d)%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 ic^db=1,2\}
            child_ig^d=2*(ig^d-1)+ic^d;
            allocate(tree%node%child(ic^d)%node)
            tree%node%child(ic^d)%node%parent%node => tree%node
            call read_node(tree%node%child(ic^d),child_ig^d,child_level)
         {end do\}
      end if
 
    end subroutine read_node
 
end subroutine read_forest