variable_halos_s.f90 Source File

submodule (variable_m) variable_halos_s
    !! Routines related with handling halos of variables 
    implicit none
contains

    module subroutine create_halos(self, comm_handler, &
                                        nhalo_fwd, nhalo_bwd)
        class(variable_t), intent(inout), target :: self  
        type(comm_handler_t), intent(in) :: comm_handler
        integer, intent(in) :: nhalo_fwd
        integer, intent(in) :: nhalo_bwd
        
        integer :: ierr, ndim_max, plane_halo, np_halo
        
        ! Number of points on all planes and maximum number of points per plane
        allocate(self%ndim_planes(0:comm_handler%get_nplanes()-1))
        call MPI_allgather(self%ndim, 1, MPI_INTEGER, &
                           self%ndim_planes, 1, MPI_INTEGER, &
                           comm_handler%get_comm_planes(), ierr)
        ndim_max = maxval(self%ndim_planes)
        
        ! Allocate communication requests
        if (allocated(self%req_fwd)) deallocate(self%req_fwd)
        if (allocated(self%req_bwd)) deallocate(self%req_bwd)
        allocate(self%req_fwd(2, nhalo_fwd))  
        allocate(self%req_bwd(2, nhalo_bwd)) 
        
        ! Allocate halo fields
        self%nhalo_fwd = nhalo_fwd    
        self%nhalo_bwd = nhalo_bwd
        if (allocated(self%halos_fwd)) deallocate(self%halos_fwd)
        if (allocated(self%halos_bwd)) deallocate(self%halos_bwd)
        allocate(self%halos_fwd(ndim_max, nhalo_fwd), source=GP_NAN)
        allocate(self%halos_bwd(ndim_max, nhalo_bwd), source=GP_NAN)
        
        ! Associate pointers to adjacent planes
        if (nhalo_fwd >= 1) then
            plane_halo = comm_handler%get_plane() + 1
            plane_halo = modulo(plane_halo, comm_handler%get_nplanes())
            np_halo = self%ndim_planes(plane_halo)
            self%hfwd => self%halos_fwd(1:np_halo,1)
        endif
        
        if (nhalo_bwd >= 1) then
            plane_halo = comm_handler%get_plane() - 1
            plane_halo = modulo(plane_halo, comm_handler%get_nplanes())
            np_halo = self%ndim_planes(plane_halo)
            self%hbwd => self%halos_bwd(1:np_halo,1)
        endif
               
    end subroutine
    
    module function get_halo(self, comm_handler, ihalo) result(res)
        real(GP), dimension(:), pointer :: res
        class(variable_t), intent(inout), target :: self 
        type(comm_handler_t), intent(in) :: comm_handler 
        integer, intent(in) :: ihalo
        
        integer :: plane_halo, np_halo
        
        ! Number of points on halo plane
        plane_halo = comm_handler%get_plane() + ihalo
        plane_halo = modulo(plane_halo, comm_handler%get_nplanes())
        np_halo = self%ndim_planes(plane_halo)
        
        if (ihalo < 0 .and. ihalo >= -self%nhalo_bwd) then
            res(1:np_halo) => self%halos_bwd(1:np_halo, -ihalo)
        elseif (ihalo > 0 .and. ihalo <= self%nhalo_fwd) then
            res(1:np_halo) => self%halos_fwd(1:np_halo, ihalo)
        elseif (ihalo == 0) then
            res(1:np_halo) => self%vals(1:np_halo)
        else
            call handle_error('ihalo not valid', &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                               error_info_t('ihalo: ',[ihalo]))
        endif
        
    end function
    
    module subroutine start_comm_halos(self, comm_handler, step)
        class(variable_t), intent(inout) :: self 
        type(comm_handler_t), intent(in) :: comm_handler
        integer, intent(in), optional :: step
    
        integer :: ndim, ndim_halo, &
                   step_start, step_end, step_current, step_abs, plane_halo
                
        ndim = self%ndim

        if (present(step)) then
            step_start = step
            step_end = step
        else 
            step_start = -self%nhalo_bwd
            step_end = self%nhalo_fwd
        endif
    
        do step_current = step_start, step_end
            step_abs = abs(step_current)
            
            plane_halo = comm_handler%get_plane() + step_current
            plane_halo = modulo(plane_halo, comm_handler%get_nplanes())
            ndim_halo = self%ndim_planes(plane_halo)
            
            if (step_current == 0) then
                cycle
            elseif (step_current < 0) then
                ! Backward communication
                call start_comm_planes(comm_handler%get_comm_planes(), &
                                    ndim, ndim_halo, self%vals, step_current, &
                                    self%halos_bwd(1:ndim_halo,step_abs), &
                                    self%req_bwd(1,step_abs) )
            elseif (step_current > 0) then
                ! Forward communication
                call start_comm_planes(comm_handler%get_comm_planes(), &
                                    ndim, ndim_halo, self%vals, step_current, &
                                    self%halos_fwd(1:ndim_halo,step_abs), &
                                    self%req_fwd(1,step_abs) )
            endif
        enddo        

    end subroutine 
    
    module subroutine finalize_comm_halos(self, comm_handler, step)
        class(variable_t), intent(inout) :: self  
        type(comm_handler_t), intent(in) :: comm_handler
        integer, intent(in), optional :: step

        integer :: step_start, step_end, step_current, step_abs

        if (present(step)) then
            step_start = step
            step_end = step
        else 
            step_start = -self%nhalo_bwd
            step_end = self%nhalo_fwd
        endif

        do step_current = step_start, step_end
            step_abs = abs(step_current)
            if (step_current == 0) then
                cycle
            elseif (step_current < 0) then
                ! Backward communication
                call finalize_comm_planes(self%req_bwd(1,step_abs))
            elseif (step_current > 0) then
                ! Forward communication
                call finalize_comm_planes(self%req_fwd(1,step_abs))
            endif
        enddo

    end subroutine  

end submodule