mod_functions_forest.t

Go to the documentation of this file.

module mod_functions_forest
 
  implicit none
  private
 
  public :: write_forest 
  public :: read_forest 
  public :: coarsen_tree_leaf 
  public :: refine_tree_leaf 
  public :: change_ipe_tree_leaf
  public :: init_forest_root
 
 
contains
 
 
  !> build root forest
  subroutine init_forest_root
    use mod_forest
    use mod_global_parameters
    use mod_space_filling_curve
    use mod_amr_solution_node, only: getnode 
    use mod_amr_neighbors, only: find_root_neighbor
 
    integer :: ig^d, level, igrid, ipe
    integer :: iside, i^d, morton_no, isfc
  
  
    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
  
    particle_mesh_epoch=particle_mesh_epoch+1_8
    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
  
    particle_mesh_epoch=particle_mesh_epoch+1_8
    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
    use mod_amr_neighbors, only: find_neighbor  
 
    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
  
    particle_mesh_epoch=particle_mesh_epoch+1_8
    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
  
    particle_mesh_epoch=particle_mesh_epoch+1_8
    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
    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 :: ig^d, level, morton_no, igrid, ipe, isfc
  
    particle_mesh_epoch=particle_mesh_epoch+1_8
    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
  
        integer :: ic^d, child_ig^d, child_level
        logical :: leaf
  
        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
 
 
end module mod_functions_forest