mod_functions_connectivity.t

Go to the documentation of this file.

module mod_functions_connectivity
 
 
  implicit none
  private
 
  public :: get_level_range
  public :: getigrids
  public :: build_connectivity
  contains
 
  subroutine get_level_range
    use mod_forest
    use mod_global_parameters
  
    integer :: level
  
    ! determine new finest level
    do level=refine_max_level,1,-1
       if (associated(level_tail(level)%node)) then
          levmax=level
          exit
       end if
    end do
  
    ! determine coarsest level
    do level=1,levmax
       if (associated(level_tail(level)%node)) then
          levmin=level
          exit
       end if
    end do
  
  end subroutine get_level_range
  
  subroutine getigrids
    use mod_forest
    use mod_global_parameters
    implicit none
  
    integer :: iigrid, igrid
  
    iigrid=0
    do igrid=1,max_blocks
       if (igrid_inuse(igrid,mype)) then
          iigrid=iigrid+1
          igrids(iigrid)=igrid
       end if
    end do
  
    igridstail=iigrid
  
  end subroutine getigrids
  
  subroutine build_connectivity
    use mod_forest
    use mod_global_parameters
    use mod_ghostcells_update
    use mod_amr_neighbors, only: find_neighbor 
  
    integer :: iigrid, igrid, i^d, my_neighbor_type
    integer :: iside, idim, ic^d, inc^d, ih^d, icdim
    type(tree_node_ptr) :: tree, my_neighbor, child
    logical, dimension(^ND) :: pole
    logical :: nopole
    ! Variables to detect special corners for stagger grid
    integer :: idir,pi^d, mi^d, ph^d, mh^d, ipe_neighbor
    integer :: nrecvs,nsends
  
    ! total size of buffer arrays
    integer :: nbuff_bc_recv_srl, nbuff_bc_send_srl, nbuff_bc_recv_r, nbuff_bc_send_r, nbuff_bc_recv_p, nbuff_bc_send_p
  
    nrecv_bc_srl=0; nsend_bc_srl=0
    nrecv_bc_r=0; nsend_bc_r=0
    nrecv_bc_p=0; nsend_bc_p=0
    nrecv_fc=0; nsend_fc=0
    nbuff_bc_recv_srl=0; nbuff_bc_send_srl=0
    nbuff_bc_recv_r=0; nbuff_bc_send_r=0
    nbuff_bc_recv_p=0; nbuff_bc_send_p=0
    if(stagger_grid) nrecv_cc=0; nsend_cc=0
  
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       tree%node => igrid_to_node(igrid,mype)%node
  
       {do i^db=-1,1\}
          ! skip the grid itself
          if (i^d==0|.and.) then
             neighbor_type(0^d&,igrid)=0
             neighbor(1,0^d&,igrid)=igrid
             neighbor(2,0^d&,igrid)=mype
          else
             call find_neighbor(my_neighbor,my_neighbor_type,tree,i^d,pole)
             nopole=.not.any(pole)
  
             select case (my_neighbor_type)
             ! adjacent to physical boundary
             case (neighbor_boundary)
                neighbor(1,i^d,igrid)=0
                neighbor(2,i^d,igrid)=-1
             ! fine-coarse transition
             case (neighbor_coarse)
                neighbor(1,i^d,igrid)=my_neighbor%node%igrid
                neighbor(2,i^d,igrid)=my_neighbor%node%ipe
                if (my_neighbor%node%ipe/=mype) then
                   ic^d=1+modulo(tree%node%ig^d-1,2);
                   if ({(i^d==0.or.i^d==2*ic^d-3)|.and.}) then
                     nrecv_bc_p=nrecv_bc_p+1
                     nsend_bc_r=nsend_bc_r+1
                     nbuff_bc_send_r=nbuff_bc_send_r+sizes_r_send_total(i^d)
                     ! This is the local index of the prolonged ghost zone
                     inc^d=ic^d+i^d;
                     nbuff_bc_recv_p=nbuff_bc_recv_p+sizes_p_recv_total(inc^d)
                   end if
                end if
             ! same refinement level
             case (neighbor_sibling)
                neighbor(1,i^d,igrid)=my_neighbor%node%igrid
                neighbor(2,i^d,igrid)=my_neighbor%node%ipe
                if (my_neighbor%node%ipe/=mype) then
                  nrecv_bc_srl=nrecv_bc_srl+1
                  nsend_bc_srl=nsend_bc_srl+1
                  nbuff_bc_send_srl=nbuff_bc_send_srl+sizes_srl_send_total(i^d)
                  nbuff_bc_recv_srl=nbuff_bc_recv_srl+sizes_srl_recv_total(i^d)
                end if
             ! coarse-fine transition
             case (neighbor_fine)
                neighbor(1,i^d,igrid)=0
                neighbor(2,i^d,igrid)=-1
                ! Loop over the local indices of children ic^D
                ! and calculate local indices of ghost zone inc^D.
                {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
                   inc^db=2*i^db+ic^db
                   if (pole(^db)) then
                      ih^db=3-ic^db
                   else
                      ih^db=ic^db
                   end if\}
                   child%node => my_neighbor%node%child(ih^d)%node
                   neighbor_child(1,inc^d,igrid)=child%node%igrid
                   neighbor_child(2,inc^d,igrid)=child%node%ipe
                   if (child%node%ipe/=mype) then
                     nrecv_bc_r=nrecv_bc_r+1
                     nsend_bc_p=nsend_bc_p+1
                     nbuff_bc_send_p=nbuff_bc_send_p+sizes_p_send_total(inc^d)
                     nbuff_bc_recv_r=nbuff_bc_recv_r+sizes_r_recv_total(inc^d)
                   end if
                {end do\}
             end select
  
             ! flux fix for conservation only for pure directional shifts
             if ({abs(i^d)+}==1) then
                {if (i^d/=0) then
                   idim=^d
                   iside=int((i^d+3)/2)
                end if\}
                select case (my_neighbor_type)
                ! only across fine-coarse or coarse-fine boundaries
                case (neighbor_coarse)
                   if (my_neighbor%node%ipe/=mype) then
                      if (.not.pole(idim)) nsend_fc(idim)=nsend_fc(idim)+1
                   end if
                case (neighbor_fine)
                   if (pole(idim)) then
                      icdim=iside
                   else
                      icdim=3-iside
                   end if
                   select case (idim)
                   {case (^d)
                      {do ic^d=icdim,icdim^d%do ic^dd=1,2\}
                         child%node => my_neighbor%node%child(ic^dd)%node
                         if (child%node%ipe/=mype) then
                            if (.not.pole(^d)) nrecv_fc(^d)=nrecv_fc(^d)+1
                         end if
                      {end do\} \}
                   end select
                end select
             end if
  
             if (phi_ > 0) then
               neighbor_pole(i^d,igrid)=0
               if (my_neighbor_type>1) then
                 do idim=1,^nd
                   if (pole(idim)) then
                     neighbor_pole(i^d,igrid)=idim
                     exit ! there can only be one pole between two meshes
                   end if
                 end do
               end if
             end if
             neighbor_type(i^d,igrid)=my_neighbor_type
  
          end if
       {end do\}
  
       if(stagger_grid) then
       !Now all the neighbour information is known.
       !Check if there are special corners that need to be communicated
       !To determine whether to send/receive, we must check three neighbours
        {do i^db=-1,1\}
           if ({abs(i^d)+}==1) then
             if (neighbor_pole(i^d,igrid)/=0) cycle
            ! Assign value to idim and iside
            {if (i^d/=0) then
                idim=^d
                iside=int((i^d+3)/2)
             end if\}
             ! Fine block surrounded by coarse blocks
             if (neighbor_type(i^d,igrid)==2) then
               do idir=idim+1,ndim
                 pi^d=i^d+kr(idir,^d);
                 mi^d=i^d-kr(idir,^d);
                 ph^d=pi^d-kr(idim,^d)*(2*iside-3);
                 mh^d=mi^d-kr(idim,^d)*(2*iside-3);
  
                 if (neighbor_type(pi^d,igrid)==2.and.&
                     neighbor_type(ph^d,igrid)==2.and.&
                     mype/=neighbor(2,pi^d,igrid).and.&
                     neighbor_pole(pi^d,igrid)==0) then
                    nsend_cc(idim) = nsend_cc(idim) + 1
                 end if
  
                 if (neighbor_type(mi^d,igrid)==2.and.&
                     neighbor_type(mh^d,igrid)==2.and.&
                     mype/=neighbor(2,mi^d,igrid).and.&
                     neighbor_pole(mi^d,igrid)==0) then
                    nsend_cc(idim) = nsend_cc(idim) + 1
                 end if
               end do
             end if
             ! Coarse block diagonal to fine block(s)
             if (neighbor_type(i^d,igrid)==3) then
               do idir=idim+1,ndim
                 pi^d=i^d+kr(idir,^d);
                 mi^d=i^d-kr(idir,^d);
                 ph^d=pi^d-kr(idim,^d)*(2*iside-3);
                 mh^d=mi^d-kr(idim,^d)*(2*iside-3);
  
                 if (neighbor_type(pi^d,igrid)==4.and.&
                     neighbor_type(ph^d,igrid)==3.and.&
                     neighbor_pole(pi^d,igrid)==0) then
                  ! Loop on children (several in 3D)
                  {do ic^db=1+int((1-pi^db)/2),2-int((1+pi^db)/2)
                      inc^db=2*pi^db+ic^db\}
                     if (mype.ne.neighbor_child(2,inc^d,igrid)) then
                       nrecv_cc(idim) = nrecv_cc(idim) + 1
                     end if
                  {end do\}
                 end if
  
                 if (neighbor_type(mi^d,igrid)==4.and.&
                     neighbor_type(mh^d,igrid)==3.and.&
                     neighbor_pole(mi^d,igrid)==0) then
                  ! Loop on children (several in 3D)
                  {do ic^db=1+int((1-mi^db)/2),2-int((1+mi^db)/2)
                      inc^db=2*mi^db+ic^db\}
                     if (mype.ne.neighbor_child(2,inc^d,igrid)) then
                       nrecv_cc(idim) = nrecv_cc(idim) + 1
                     end if
                  {end do\}
                 end if
               end do
             end if
           end if
        {end do\}
       end if
  
    end do
  
    ! allocate space for mpi recieve for siblings and restrict ghost cell filling
    nrecvs=nrecv_bc_srl+nrecv_bc_r
    if (allocated(recvstatus_c_sr)) then
      deallocate(recvstatus_c_sr,recvrequest_c_sr)
      allocate(recvstatus_c_sr(mpi_status_size,nrecvs),recvrequest_c_sr(nrecvs))
    else
      allocate(recvstatus_c_sr(mpi_status_size,nrecvs),recvrequest_c_sr(nrecvs))
    end if
    recvrequest_c_sr=mpi_request_null
  
    ! allocate space for mpi send for siblings and restrict ghost cell filling
    nsends=nsend_bc_srl+nsend_bc_r
    if (allocated(sendstatus_c_sr)) then
      deallocate(sendstatus_c_sr,sendrequest_c_sr)
      allocate(sendstatus_c_sr(mpi_status_size,nsends),sendrequest_c_sr(nsends))
    else
      allocate(sendstatus_c_sr(mpi_status_size,nsends),sendrequest_c_sr(nsends))
    end if
    sendrequest_c_sr=mpi_request_null
  
    ! allocate space for mpi recieve for prolongation ghost cell filling
    if (allocated(recvstatus_c_p)) then
      deallocate(recvstatus_c_p,recvrequest_c_p)
      allocate(recvstatus_c_p(mpi_status_size,nrecv_bc_p),recvrequest_c_p(nrecv_bc_p))
    else
      allocate(recvstatus_c_p(mpi_status_size,nrecv_bc_p),recvrequest_c_p(nrecv_bc_p))
    end if
    recvrequest_c_p=mpi_request_null
  
    ! allocate space for mpi send for prolongation ghost cell filling
    if (allocated(sendstatus_c_p)) then
      deallocate(sendstatus_c_p,sendrequest_c_p)
      allocate(sendstatus_c_p(mpi_status_size,nsend_bc_p),sendrequest_c_p(nsend_bc_p))
    else
      allocate(sendstatus_c_p(mpi_status_size,nsend_bc_p),sendrequest_c_p(nsend_bc_p))
    end if
    sendrequest_c_p=mpi_request_null
  
    if(stagger_grid) then
      ! allocate space for recieve buffer for siblings ghost cell filling
      if (allocated(recvbuffer_srl)) then
        if (nbuff_bc_recv_srl /= size(recvbuffer_srl)) then
          deallocate(recvbuffer_srl)
          allocate(recvbuffer_srl(nbuff_bc_recv_srl))
        end if
      else
        allocate(recvbuffer_srl(nbuff_bc_recv_srl))
      end if
      if (allocated(recvstatus_srl)) then
        deallocate(recvstatus_srl,recvrequest_srl)
        allocate(recvstatus_srl(mpi_status_size,nrecv_bc_srl),recvrequest_srl(nrecv_bc_srl))
      else
        allocate(recvstatus_srl(mpi_status_size,nrecv_bc_srl),recvrequest_srl(nrecv_bc_srl))
      end if
      recvrequest_srl=mpi_request_null
  
      ! allocate space for send buffer for siblings ghost cell filling
      if (allocated(sendbuffer_srl)) then
        if (nbuff_bc_send_srl /= size(sendbuffer_srl)) then
          deallocate(sendbuffer_srl)
          allocate(sendbuffer_srl(nbuff_bc_send_srl))
        end if
      else
        allocate(sendbuffer_srl(nbuff_bc_send_srl))
      end if
      if (allocated(sendstatus_srl)) then
        deallocate(sendstatus_srl,sendrequest_srl)
        allocate(sendstatus_srl(mpi_status_size,nsend_bc_srl),sendrequest_srl(nsend_bc_srl))
      else
        allocate(sendstatus_srl(mpi_status_size,nsend_bc_srl),sendrequest_srl(nsend_bc_srl))
      end if
      sendrequest_srl=mpi_request_null
  
      ! allocate space for recieve buffer for restrict ghost cell filling
      if (allocated(recvbuffer_r)) then
        if (nbuff_bc_recv_r /= size(recvbuffer_r)) then
          deallocate(recvbuffer_r)
          allocate(recvbuffer_r(nbuff_bc_recv_r))
        end if
      else
        allocate(recvbuffer_r(nbuff_bc_recv_r))
      end if
      if (allocated(recvstatus_r)) then
        deallocate(recvstatus_r,recvrequest_r)
        allocate(recvstatus_r(mpi_status_size,nrecv_bc_r),recvrequest_r(nrecv_bc_r))
      else
        allocate(recvstatus_r(mpi_status_size,nrecv_bc_r),recvrequest_r(nrecv_bc_r))
      end if
      recvrequest_r=mpi_request_null
  
      ! allocate space for send buffer for restrict ghost cell filling
      if (allocated(sendbuffer_r)) then
        if (nbuff_bc_send_r /= size(sendbuffer_r)) then
          deallocate(sendbuffer_r)
          allocate(sendbuffer_r(nbuff_bc_send_r))
        end if
      else
        allocate(sendbuffer_r(nbuff_bc_send_r))
      end if
      if (allocated(sendstatus_r)) then
        deallocate(sendstatus_r,sendrequest_r)
        allocate(sendstatus_r(mpi_status_size,nsend_bc_r),sendrequest_r(nsend_bc_r))
      else
        allocate(sendstatus_r(mpi_status_size,nsend_bc_r),sendrequest_r(nsend_bc_r))
      end if
      sendrequest_r=mpi_request_null
  
      ! allocate space for recieve buffer for prolong ghost cell filling
      if (allocated(recvbuffer_p)) then
        if (nbuff_bc_recv_p /= size(recvbuffer_p)) then
          deallocate(recvbuffer_p)
          allocate(recvbuffer_p(nbuff_bc_recv_p))
        end if
      else
        allocate(recvbuffer_p(nbuff_bc_recv_p))
      end if
      if (allocated(recvstatus_p)) then
        deallocate(recvstatus_p,recvrequest_p)
        allocate(recvstatus_p(mpi_status_size,nrecv_bc_p),recvrequest_p(nrecv_bc_p))
      else
        allocate(recvstatus_p(mpi_status_size,nrecv_bc_p),recvrequest_p(nrecv_bc_p))
      end if
      recvrequest_p=mpi_request_null
  
      ! allocate space for send buffer for restrict ghost cell filling
      if (allocated(sendbuffer_p)) then
        if (nbuff_bc_send_p /= size(sendbuffer_p)) then
          deallocate(sendbuffer_p)
          allocate(sendbuffer_p(nbuff_bc_send_p))
        end if
      else
        allocate(sendbuffer_p(nbuff_bc_send_p))
      end if
      if (allocated(sendstatus_p)) then
        deallocate(sendstatus_p,sendrequest_p)
        allocate(sendstatus_p(mpi_status_size,nsend_bc_p),sendrequest_p(nsend_bc_p))
      else
        allocate(sendstatus_p(mpi_status_size,nsend_bc_p),sendrequest_p(nsend_bc_p))
      end if
      sendrequest_p=mpi_request_null
    end if
  
  end subroutine build_connectivity
 
end module mod_functions_connectivity