psmile_get_halo_indices.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PSMILe_get_halo_indices
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_get_halo_indices ( comp_id, grid_id_list, ierror )
00012 !
00013 ! !USES:
00014 !
00015       use PRISM_constants
00016 !
00017       use PSMILe, dummy_interface => PSMILe_get_halo_indices
00018 
00019       implicit none
00020 !
00021 ! !INPUT PARAMETERS:
00022 !
00023       Integer, Intent (In)         :: comp_id
00024 
00025       Integer, Intent (In)         :: grid_id_list(Number_of_Grids_allocated)
00026 
00027 !
00028 ! !OUTPUT PARAMETERS:
00029 !
00030       Integer, Intent (Out)        :: ierror
00031 
00032 !     Returns the error code of prism_enddef;
00033 !             ierror = 0 : No error
00034 !             ierror > 0 : Severe error
00035 !
00036 ! !LOCAL VARIABLES
00037 !
00038       Integer                      :: i, j, n
00039       Integer                      :: ii, jj
00040       Integer                      :: nbr_blocks
00041       Integer, Parameter           :: nbr_halos_inc = 8
00042       Integer                      :: halo_index
00043       Integer                      :: len
00044       Integer                      :: npoints
00045       Integer                      :: grid_id
00046       Integer                      :: comm
00047       Integer                      :: rank
00048       Integer                      :: npes
00049       Integer                      :: all_parts
00050       Integer                      :: ngrids
00051       Integer                      :: halo_range(2,ndim_3d)
00052       Integer                      :: halo_send_range(2,ndim_3d)
00053 
00054       Integer, Allocatable         :: all_ngrids(:)
00055       Integer, Allocatable         :: grange(:,:,:), all_grange(:,:,:)
00056       Integer, Allocatable         :: displs(:)
00057       Integer, Allocatable         :: counts(:)
00058 
00059       Type (Halo_info), Pointer    :: halo_info_p
00060       Type (Halo_info), Pointer    :: halo (:)
00061       Type (Halo_info), Pointer    :: new_halo (:)
00062 
00063       Integer, Parameter           :: nerrp = 2
00064       Integer                      :: ierrp (nerrp)
00065 
00066 #ifdef DEBUG
00067       Integer                      :: grid_id_old
00068 #endif
00069 !
00070 ! !DESCRIPTION:
00071 !
00072 ! Subroutine "PSMILe_get_halo_indices" exchanges halo points of with
00073 ! neighbour processes. Relevent on source processes only.
00074 !
00075 ! TODO:
00076 !
00077 ! If there is more than one component per process we have to proceed in
00078 ! two steps. Frist collect the data, then evaluate.
00079 !
00080 ! !REVISION HISTORY:
00081 !   Date      Programmer   Description
00082 ! ----------  ----------   -----------
00083 ! 08.12.01    R. Redler    created
00084 !
00085 !EOP
00086 !----------------------------------------------------------------------
00087 !
00088 ! $Id: psmile_get_halo_indices.F90 2785 2010-11-29 15:15:19Z redler $
00089 ! $Author: redler $
00090 !
00091    Character(len=len_cvs_string), save :: mycvs = 
00092        '$Id: psmile_get_halo_indices.F90 2785 2010-11-29 15:15:19Z redler $'
00093 !
00094 !----------------------------------------------------------------------
00095 !
00096 !  Initialization
00097 !
00098 #ifdef VERBOSE
00099       print 9990, trim(ch_id)
00100 
00101       call psmile_flushstd
00102 #endif /* VERBOSE */
00103 
00104       ierror  = 0
00105 
00106       comm    = Comps(comp_id)%comm_halo
00107 
00108       ! we still use the rank in the component communicator
00109 
00110       rank    = Comps(comp_id)%rank
00111 
00112       call MPI_Comm_size (Comps(comp_id)%comm_halo, npes, ierror)
00113 
00114       ! npes    = Comps(comp_id)%size
00115 
00116       halo_index = -1
00117       ngrids = 0
00118 
00119       do n = 1, Number_of_Grids_allocated
00120          if ( grid_id_list(n) /= PSMILe_Undef ) then
00121              nbr_blocks = size(Grids(n)%partition(:,1))
00122              ngrids = ngrids + nbr_blocks
00123          endif
00124       enddo
00125 
00126 #ifdef DEBUG
00127       PRINT*, 'Number of grids with blocks :', ngrids
00128       call psmile_flushstd
00129 #endif
00130 
00131       Allocate(all_ngrids(npes), STAT=ierror)
00132       if (ierror /= 0) then
00133          ierrp (1) = ierror
00134          ierrp (2) = npes
00135          ierror = PRISM_Error_Alloc
00136          call psmile_error ( ierror, 'all_ngrids', ierrp, 2, &
00137               __FILE__, __LINE__ )
00138          return
00139       endif
00140 
00141       all_parts = 0
00142 
00143       call MPI_Allgather ( ngrids, 1, MPI_Integer, &
00144            all_ngrids, 1, MPI_Integer, comm, ierror )
00145 
00146       all_parts = sum(all_ngrids)
00147 
00148       if ( all_parts == 0 ) goto 100
00149 
00150       Allocate(grange(2,ndim_3d+2, ngrids), STAT=ierror)
00151       if (ierror /= 0) then
00152          ierrp (1) = ierror
00153          ierrp (2) = len
00154          ierror = PRISM_Error_Alloc
00155          call psmile_error ( ierror, 'grange', ierrp, 2, &
00156               __FILE__, __LINE__ )
00157          return
00158       endif
00159       !
00160       ! Loop over grids and partitions to collect global information
00161       !
00162       ngrids = 0
00163 
00164       do grid_id = 1, Number_of_Grids_allocated
00165 
00166          ! To ensure that we work for the specified component only
00167 
00168          if ( Grids(grid_id)%comp_id /= comp_id ) cycle
00169 
00170          if ( grid_id_list(grid_id) /= PSMILe_Undef ) then
00171 #ifdef PRISM_ASSERTION
00172             if ( .not. associated(Grids(grid_id)%partition) ) then
00173                call psmile_assert (__FILE__, __LINE__, &
00174                     'Routine requires user call to prism_def_partition.')
00175             endif
00176 #endif
00177             nbr_blocks = size(Grids(grid_id)%partition(:,1))
00178             do n = 1, nbr_blocks
00179                ngrids = ngrids + 1
00180                do i = 1, ndim_3d
00181                   grange(1,i,ngrids) = Grids(grid_id)%partition(n,i) + 1
00182                   grange(2,i,ngrids) = Grids(grid_id)%partition(n,i) + Grids(grid_id)%extent(n,i)
00183                enddo
00184                grange(1,ndim_3d+1,ngrids) = grid_id
00185                grange(2,ndim_3d+1,ngrids) = rank
00186                grange(1,ndim_3d+2,ngrids) = n
00187                grange(2,ndim_3d+2,ngrids) = Grids(grid_id)%global_grid_id
00188             enddo
00189          endif
00190       enddo
00191 
00192       Allocate(displs(npes), counts(npes), STAT=ierror)
00193       if (ierror /= 0) then
00194          ierrp (1) = ierror
00195          ierrp (2) = 2*npes
00196          ierror = PRISM_Error_Alloc
00197          call psmile_error ( ierror, 'displs and counts', ierrp, 2, &
00198               __FILE__, __LINE__ )
00199          return
00200       endif
00201 
00202       Allocate(all_grange(2, ndim_3d+2, all_parts), STAT=ierror)
00203       if (ierror /= 0) then
00204          ierrp (1) = ierror
00205          ierrp (2) = 2 * ( ndim_3d+2 ) * all_parts
00206          ierror = PRISM_Error_Alloc
00207          call psmile_error ( ierror, 'all_grange', ierrp, 2, &
00208               __FILE__, __LINE__ )
00209          return
00210       endif
00211 
00212       displs = 0
00213       counts = 0
00214 
00215       len = 2*(ndim_3d+2)
00216 
00217       displs(1) = 0
00218       do i = 2, npes
00219          displs(i) = displs(i-1)+all_ngrids(i-1)*len
00220       enddo
00221 
00222       counts(1:npes) = all_ngrids(1:npes)*len
00223 
00224       call MPI_Allgatherv (grange, len*ngrids, MPI_Integer, &
00225            all_grange, counts, displs, MPI_Integer, &
00226            comm, ierror )
00227 
00228       Deallocate(displs, counts, STAT=ierror)
00229       if (ierror /= 0) then
00230          ierrp (1) = ierror
00231          ierrp (2) = 2*npes
00232          ierror = PRISM_Error_Dealloc
00233          call psmile_error ( ierror, 'displs and counts', ierrp, 2, &
00234               __FILE__, __LINE__ )
00235          return
00236       endif
00237 
00238 #ifdef DEBUG
00239       print *, "gathered grid information:"
00240       print *, "Grid #|global index range                  |grid id|   rank|block #|global grid id"
00241       do i = 1, all_parts
00242          print '(i7,"|",6i6,"|",3(i7,"|"),i8)', i, all_grange(:,:,i)
00243       enddo
00244 #endif
00245       !
00246       ! Determine where to send and from where to receive data
00247       ! The loop over all_part can become quite expensive, but
00248       ! as we currently lack any information about the connectivity
00249       ! of processes we do not have any other chance for the moment.
00250       !
00251 
00252       do n = 1, ngrids
00253 
00254          grid_id = grange(1,ndim_3d+1,n)
00255 
00256          if ( .not. associated(Grids(grid_id)%halo) )  then
00257             Allocate(Grids(grid_id)%halo(nbr_halos_inc), STAT=ierror)
00258             if (ierror /= 0) then
00259                ierrp (1) = ierror
00260                ierrp (2) = nbr_halos_inc
00261                ierror = PRISM_Error_Alloc
00262                call psmile_error ( ierror, 'halo', ierrp, 2, &
00263                     __FILE__, __LINE__ )
00264                return
00265             endif
00266 
00267          endif
00268 
00269          select case (Grids(grid_id)%grid_type)
00270 
00271          case ( PRISM_Irrlonlatvrt )
00272 
00273 #ifdef TODO
00274             do i = 1, ndim_3d
00275 
00276                ! determine ranks for halo exchange
00277 
00278                if ( (grange(1,i,n)-1 >= all_grange(1,i,j)  .and. &
00279                      grange(1,i,n)-1 <= all_grange(2,i,j)) .or.  &
00280                     (grange(2,i,n)+1 >= all_grange(1,i,j)  .and. &
00281                      grange(2,i,n)+1 <= all_grange(2,i,j)) ) then
00282 
00283                   ii = mod(i  ,ndim_3d)+1
00284                   jj = mod(i+1,ndim_3d)+1
00285 
00286                   npoints = ( min(grange(2,ii,n),all_grange(2,ii,j)) - &
00287                               max(grange(1,ii,n),all_grange(1,ii,j)) + 1 ) * &
00288                             ( min(grange(2,jj,n),all_grange(2,jj,j)) - &
00289                               max(grange(1,jj,n),all_grange(1,jj,j)) + 1 )
00290 
00291                   print * , ' exchange ', npoints, ' with rank ', all_grange(2,ndim_3d+1,j)
00292 
00293                endif
00294             enddo
00295 #endif
00296             ierrp (1) = grid_id
00297             ierrp (2) = Grids(grid_id)%grid_type
00298 
00299             ierror = PRISM_Error_Internal
00300 
00301             call psmile_error ( ierror, 'grid type not yet supported for halo exchange', &
00302                  ierrp, 2, __FILE__, __LINE__ )
00303 
00304          case ( PRISM_Irrlonlat_regvrt )
00305 
00306             do j = 1, all_parts
00307 
00308                if ( Grids(grid_id)%global_grid_id /= grange(2,ndim_3d+2,n) ) cycle
00309 
00310                ! if the currently range is on a remote processes
00311                ! MoHa: what about another block of the same grid on
00312                ! the local process, shouldn't we check it as well?
00313                if ( all_grange(2,ndim_3d+1,j) /= rank ) then
00314 
00315                   do i = 1, ndim_2d
00316 
00317                      ! determine ranks for halo exchange
00318 
00319                      ! if the local halo is within the index range the the
00320                      ! currently tested grid range
00321                      if ( (grange(1,i,n)-1 >= all_grange(1,i,j)  .and. &
00322                            grange(1,i,n)-1 <= all_grange(2,i,j)) .or.  &
00323                           (grange(2,i,n)+1 >= all_grange(1,i,j)  .and. &
00324                            grange(2,i,n)+1 <= all_grange(2,i,j)) ) then
00325 
00326                         ! check whether the grid ranges match in the other grid dimension
00327                         ii = mod(i  ,ndim_2d)+1
00328                         if ( grange(2,ii,n) < all_grange(1,ii,j) .or. &
00329                              grange(1,ii,n) > all_grange(2,ii,j) ) cycle
00330 
00331                         Grids(grid_id)%nbr_halo_segments = &
00332                            Grids(grid_id)%nbr_halo_segments + 1
00333 
00334                         halo_index = Grids(grid_id)%nbr_halo_segments
00335 
00336                         ! check which halo part matches with the tested grid range
00337                         ! the halo we need locally is halo_range
00338                         ! once we have halo_range we automatically know, that the
00339                         ! other process needs matching halo from the local data (halo_send_range)
00340                         if ( grange(1,i,n)-1 >= all_grange(1,i,j) .and. &
00341                              grange(1,i,n)-1 <= all_grange(2,i,j) ) then
00342 
00343                            ! ... lower bound
00344 
00345                            if ( i == 1 ) then
00346 
00347                                halo_range(1:2,1) = grange(1,1,n)-1
00348                                halo_range(1,2) = max(grange(1,2,n),all_grange(1,2,j))
00349                                halo_range(2,2) = min(grange(2,2,n),all_grange(2,2,j))
00350 
00351                                halo_send_range(1:2,1) = halo_range(1:2,1) + 1
00352                                halo_send_range(1:2,2) = halo_range(1:2,2)
00353 
00354                            else
00355 
00356                                halo_range(1,1) = max(grange(1,1,n),all_grange(1,1,j))
00357                                halo_range(2,1) = min(grange(2,1,n),all_grange(2,1,j))
00358                                halo_range(1:2,2) = grange(1,2,n)-1
00359 
00360                                halo_send_range(1:2,1) = halo_range(1:2,1)
00361                                halo_send_range(1:2,2) = halo_range(1:2,2) + 1
00362 
00363                            endif
00364 
00365                         else
00366 
00367                            ! ... upper bound
00368 
00369                            if ( i == 1 ) then
00370 
00371                                halo_range(1:2,1) = grange(2,1,n)+1
00372                                halo_range(1,2) = max(grange(1,2,n),all_grange(1,2,j))
00373                                halo_range(2,2) = min(grange(2,2,n),all_grange(2,2,j))
00374 
00375                                halo_send_range(1:2,1) = halo_range(1:2,1)-1
00376                                halo_send_range(1:2,2) = halo_range(1:2,2)
00377 
00378                            else
00379 
00380                                halo_range(1,1) = max(grange(1,1,n),all_grange(1,1,j))
00381                                halo_range(2,1) = min(grange(2,1,n),all_grange(2,1,j))
00382                                halo_range(1:2,2) = grange(2,2,n)+1
00383 
00384                                halo_send_range(1:2,1) = halo_range(1:2,1)
00385                                halo_send_range(1:2,2) = halo_range(1:2,2)-1
00386 
00387                            endif
00388 
00389                         endif
00390 
00391                         npoints = (halo_range(2,2)-halo_range(1,2)+1) * &
00392                                   (halo_range(2,1)-halo_range(1,1)+1)
00393 
00394                         halo_info_p => Grids(grid_id)%halo(halo_index)
00395 
00396                         halo_info_p%local_rank            = rank
00397                         halo_info_p%remote_rank           = all_grange(2,ndim_3d+1,j)
00398                         halo_info_p%remote_grid_id        = all_grange(1,ndim_3d+1,j)
00399                         halo_info_p%remote_block          = all_grange(1,ndim_3d+2,j)
00400                         halo_info_p%halo_size             = npoints
00401                         halo_info_p%local_block           = grange(1,ndim_3d+2,n)
00402                         halo_info_p%global_range(1:2,1:2) = halo_range(1:2,1:2)
00403                         halo_info_p%global_range(1:2,3)   = 1
00404                         halo_info_p%send_range(1:2,1:2)   = halo_send_range(1:2,1:2)
00405                         halo_info_p%send_range(1:2,3)     = 1
00406 
00407 #ifdef DEBUG
00408                         print *, "computed halo exchange information:"
00409                         print *, "(l == local, r == remote)"
00410                         WRITE(*,7770)
00411                         WRITE(*,7771) halo_info_p%local_rank, &
00412                                    halo_info_p%remote_rank, &
00413                                    grid_id, &
00414                                    halo_info_p%remote_grid_id, &
00415                                    halo_info_p%local_block, &
00416                                    halo_info_p%remote_block, &
00417                                    halo_info_p%halo_size, &
00418                                    halo_info_p%global_range
00419 
00420                         WRITE(*,7772)
00421                         WRITE(*,7771) halo_info_p%local_rank, &
00422                                     halo_info_p%remote_rank, &
00423                                     grid_id, &
00424                                     halo_info_p%remote_grid_id, &
00425                                     halo_info_p%local_block, &
00426                                     halo_info_p%remote_block, &
00427                                     halo_info_p%halo_size, &
00428                                     halo_info_p%send_range
00429 #endif
00430                      endif
00431 
00432                   enddo
00433 
00434                endif
00435 
00436             enddo
00437 
00438          case DEFAULT
00439 
00440             ierrp (1) = grid_id
00441             ierrp (2) = Grids(grid_id)%grid_type
00442 
00443             ierror = PRISM_Error_Internal
00444 
00445             call psmile_error ( ierror, 'grid type not yet supported for halo exchange', &
00446                  ierrp, 2, __FILE__, __LINE__ )
00447 
00448          end select
00449 
00450      enddo ! n
00451 #ifdef DEBUG
00452       DO n = 1, ngrids
00453          grid_id = grange(1,ndim_3d+1,n)
00454          IF (n == 1) THEN
00455              PRINT*, 'In get_halo_indices nbr_halo_segments :',Grids(grid_id)%nbr_halo_segments,&
00456                 'for grid_id :',grid_id
00457          ELSE
00458              IF (grid_id .NE. grid_id_old) THEN
00459                  PRINT*, 'In get_halo_indices nbr_halo_segments :',Grids(grid_id)%nbr_halo_segments,&
00460                     'for grid_id :',grid_id
00461              ENDIF
00462          ENDIF
00463         grid_id_old=grid_id
00464        ENDDO
00465 #endif
00466 
00467 
00468       ! The loop over n works on all grid ids and partitions. This mix gets sorted
00469       ! into grid ids. Now we shrink the Grids(grid_id)%halo to a size equal to
00470       ! halo_index for each gid id.
00471 
00472       do grid_id = 1, Number_of_Grids_allocated
00473 
00474          if ( associated( Grids(grid_id)%halo ) ) then
00475 
00476             halo => Grids(grid_id)%halo
00477             halo_index = Grids(grid_id)%nbr_halo_segments
00478 
00479             if (halo_index > 0 .and. halo_index /= size(halo)) then
00480 
00481                Allocate(new_halo(halo_index), STAT=ierror)
00482                if (ierror /= 0) then
00483                   ierrp (1) = ierror
00484                   ierrp (2) = halo_index
00485                   ierror = PRISM_Error_Alloc
00486                   call psmile_error ( ierror, 'halo', ierrp, 2, &
00487                         __FILE__, __LINE__ )
00488                   return
00489                endif
00490 
00491                new_halo(:) = halo(1:halo_index)
00492 
00493                Deallocate(halo, STAT=ierror)
00494                if (ierror /= 0) then
00495                   ierrp (1) = ierror
00496                   ierrp (2) = nbr_halos_inc
00497                   ierror = PRISM_Error_Dealloc
00498                   call psmile_error ( ierror, 'halo', ierrp, 2, &
00499                         __FILE__, __LINE__ )
00500                   return
00501                endif
00502 
00503                Grids(grid_id)%halo => new_halo
00504                Nullify(new_halo)
00505 
00506             endif
00507          endif
00508       enddo
00509 
00510       !
00511       ! -----------------------------------------------------------------
00512       ! Clean up memory
00513       ! -----------------------------------------------------------------
00514       !
00515 
00516       Deallocate(all_grange, STAT=ierror)
00517       if (ierror /= 0) then
00518          ierrp (1) = ierror
00519          ierrp (2) = 2 * ( ndim_3d+2 ) * all_parts
00520          ierror = PRISM_Error_Dealloc
00521          call psmile_error ( ierror, 'all_grange', ierrp, 2, &
00522                __FILE__, __LINE__ )
00523          return
00524       endif
00525 
00526 100  continue
00527 
00528       Deallocate(all_ngrids, STAT=ierror)
00529       if (ierror /= 0) then
00530          ierrp (1) = ierror
00531          ierrp (2) = npes
00532          ierror = PRISM_Error_Dealloc
00533          call psmile_error ( ierror, 'all_ngrids', ierrp, 2, &
00534                __FILE__, __LINE__ )
00535          return
00536       endif
00537 !
00538 #ifdef VERBOSE
00539       print 9980, trim(ch_id), ierror
00540 
00541       call psmile_flushstd
00542 #endif /* VERBOSE */
00543 !
00544 !  Formats:
00545 !
00546 7770 format ("l rank|r rank|l grid id|r grid id|l block #|r block #|halo size|range recv")
00547 7771 format (2(i6,"|"),5(i9,"|"),6i6)
00548 7772 format ("l rank|r rank|l grid id|r grid id|l block #|r block #|halo size|range sent")
00549 9990 format (1x, a, ': psmile_get_halo_indices')
00550 9980 format (1x, a, ': psmile_get_halo_indices eof: ierror ', i3)
00551 
00552       end subroutine PSMILe_get_halo_indices
00553