psmile_get_halo_points.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! Copyright 2010, DKRZ, Hamburg, Germany.
00004 ! All rights reserved. Use is subject to OASIS4 license terms.
00005 !-----------------------------------------------------------------------
00006 !BOP
00007 !
00008 ! !ROUTINE: psmile_get_halo_points
00009 !
00010 ! !INTERFACE:
00011 
00012 subroutine psmile_get_halo_points (comp_id, ierror)
00013 !
00014 ! !USES:
00015 !
00016    use psmile, only            : halo_info, halo_block, fields, grids, methods, comps, &
00017                                  coords_block
00018 
00019    use psmile_common, only     : ndim_2d, ndim_3d, ch_id, number_of_methods_allocated,     &
00020                                  psmile_status_defined, psmile_undef, len_cvs_string,      &
00021                                  MPI_Request_Null, MPI_STATUS_SIZE, MPI_REAL, &
00022                                  MPI_DOUBLE_PRECISION, halotag
00023 
00024    use prism_constants, only   : prism_irrlonlat_regvrt, &
00025                                  prism_irrlonlatvrt,     &
00026                                  prism_error_alloc
00027 
00028    use psmile_reallocate, only : psmile_realloc
00029 
00030    implicit none
00031 !
00032 ! !INPUT PARAMETERS:
00033 !
00034    integer, intent (in)         :: comp_id
00035 !
00036 ! !OUTPUT PARAMETERS:
00037 !
00038    integer, intent (out)        :: ierror
00039 
00040 !     Returns the error code of prism_enddef;
00041 !             ierror = 0 : No error
00042 !             ierror > 0 : Severe error
00043 !
00044 ! !LOCAL VARIABLES
00045 !
00046    integer                   :: grid_id
00047    integer                   :: halo_id
00048    integer                   :: field_id
00049    integer                   :: method_id
00050    integer                   :: shift(ndim_3d) ! required to transform global indices into local ones
00051    type (halo_info), pointer :: halo
00052 
00053    logical                   :: method_needs_halos(number_of_methods_allocated)
00054    integer                   :: num_requests
00055    integer, pointer          :: requests(:)
00056    integer, allocatable      :: statuses(:,:)
00057    double precision, pointer :: dble_send_buf(:)
00058    real, pointer             :: real_send_buf(:)
00059    integer                   :: ndim
00060 
00061 !
00062 ! !DESCRIPTION:
00063 !
00064 ! Subroutine "psmile_get_halo_points" exchanges halo points of with
00065 ! neighbour processes. Relevent on source processes only.
00066 !
00067 !
00068 ! !REVISION HISTORY:
00069 !   Date      Programmer   Description
00070 ! ----------  ----------   -----------
00071 ! 08.12.01    R. Redler    created
00072 ! 16.11.10    M. Hanke     rewrote psmile_get_halo_points based on previouse version
00073 !
00074 !EOP
00075 !----------------------------------------------------------------------
00076 !
00077 ! $Id: psmile_get_halo_points.F90 2790 2010-12-01 11:26:54Z coquart $
00078 ! $Author: coquart $
00079 !
00080    Character(len=len_cvs_string), save :: mycvs = 
00081        '$Id: psmile_get_halo_points.F90 2790 2010-12-01 11:26:54Z coquart $'
00082 !
00083 ! ---------------------------------------------------------------------
00084 !
00085 !  initialization
00086 !
00087 #ifdef VERBOSE
00088    print 9990, trim(ch_id)
00089    call psmile_flushstd
00090 #endif /* VERBOSE */
00091 
00092    ierror  = 0
00093    num_requests = 0
00094 
00095    nullify (requests)
00096    nullify (dble_send_buf)
00097    nullify (real_send_buf)
00098 
00099    ! for all grids
00100    do grid_id = 1, size (Grids)
00101 
00102       if ( Grids(grid_id)%status /= psmile_status_defined .or.      &
00103            Grids(grid_id)%comp_id /= comp_id .or.                   &
00104           (Grids(grid_id)%grid_type /= prism_irrlonlat_regvrt .and. &
00105            Grids(grid_id)%grid_type /= prism_irrlonlatvrt) .or.     &
00106            Grids(grid_id)%nbr_halo_segments < 1 ) then
00107          cycle
00108       endif
00109 
00110 #ifdef PRISM_ASSERTION
00111       if ( .not. associated (Grids(grid_id)%partition)) then
00112          call psmile_assert (__FILE__, __LINE__, &
00113                              "psmile_get_halo_points requires partition information")
00114       endif
00115       if (size (Grids(grid_id)%partition, 1) /= 1) then
00116          call psmile_assert (__FILE__, __LINE__, &
00117                              "psmile_get_halo_points does not support more than one block per grid")
00118       endif
00119 #endif
00120 
00121       method_needs_halos = .false.
00122 
00123       ! get the methods for which the halos need to be exchanged
00124 
00125       do field_id = 1, size (Fields)
00126 
00127          ! if we have a valid field, if it is a source
00128 
00129          if (fields(field_id)%smioc_loc /= psmile_undef .and. &
00130              associated (fields(field_id)%taskout)) then
00131 
00132             ! if the grid ids match
00133 
00134             if ( methods(fields(field_id)%method_id)%grid_id == grid_id) then
00135 
00136                method_needs_halos(fields(field_id)%method_id) = .true.
00137 
00138                call init_method_halos(fields(field_id)%method_id, &
00139                                       Grids(grid_id)%nbr_halo_segments)
00140             endif
00141          endif
00142       enddo ! field_id
00143 
00144       ! if there are no methods we have to work on
00145 
00146       if (count (method_needs_halos) == 0) cycle
00147 
00148       ! compute offset required to transform the stored halo indices
00149       ! (in global coordinates) into the local index space
00150 
00151       shift = Grids(grid_id)%grid_shape(1,:)-(Grids(grid_id)%partition(1,:) + 1)
00152 
00153       ! get the number of dimension which are going to be exchanged
00154 
00155       if ( Grids(grid_id)%grid_type == prism_irrlonlat_regvrt ) then
00156          ndim = ndim_2d
00157       else
00158          ndim = ndim_3d
00159       endif
00160 
00161       ! allocate request objects for receiving the halos
00162 
00163       requests => psmile_realloc(requests, num_requests + Grids(grid_id)%nbr_halo_segments &
00164                                  * count (method_needs_halos) * ndim)
00165       requests(num_requests+1:) = MPI_Request_Null
00166 
00167       ! for all halo segments of the current grid
00168 
00169       do halo_id = 1, Grids(grid_id)%nbr_halo_segments
00170 
00171          ! get halo
00172          halo => Grids(grid_id)%halo(halo_id)
00173 
00174          ! for all methods associated with current grid
00175 
00176          do method_id = 1, size (method_needs_halos)
00177 
00178             if (method_needs_halos(method_id)) then
00179 
00180                ! there is always a matching halo on the other process,
00181                ! therefore we can send and request one at the same time
00182 
00183                ! post request for halo segment
00184 
00185                call recv_halo (method_id, halo, halo_id, shift, ndim, &
00186                               requests(num_requests+1:num_requests+ndim), ierror)
00187                num_requests = num_requests + ndim
00188 
00189                ! bsend halo to associated other process
00190 
00191                call send_halo (method_id, halo, shift, ndim, ierror)
00192 
00193             endif ! method_needs_halos(i)
00194          enddo ! method_id
00195       enddo ! halo_id
00196    enddo ! grid_id
00197 
00198    if (num_requests > 0) then
00199 
00200       allocate (statuses(MPI_STATUS_SIZE,num_requests), stat=ierror)
00201       if (ierror /= 0) then
00202          ierror = prism_error_alloc
00203          call psmile_error (ierror, 'statuses', (/num_requests * MPI_STATUS_SIZE/), &
00204                            1, __FILE__, __LINE__)
00205          return
00206       endif
00207 
00208       ! waitall
00209 
00210       call MPI_Waitall ( num_requests, requests, statuses, ierror )
00211 
00212       deallocate (requests, statuses)
00213 
00214       if (associated (dble_send_buf)) deallocate (dble_send_buf)
00215       if (associated (real_send_buf)) deallocate (real_send_buf)
00216 
00217    endif
00218 
00219 #ifdef VERBOSE
00220    print 9980, trim(ch_id), ierror
00221    call psmile_flushstd
00222 #endif /* VERBOSE */
00223 !
00224 !  Formats:
00225 !
00226 9990 format (1x, a, ': psmile_get_halo_points')
00227 9980 format (1x, a, ': psmile_get_halo_points; eof ierror =', i4)
00228 
00229 contains
00230 
00231    ! ===========================================================================
00232 
00233    subroutine init_method_halos (method_id, num_halos)
00234 
00235       integer, intent (in) :: method_id
00236       integer, intent (in) :: num_halos
00237 
00238       integer :: i
00239       integer :: j
00240 
00241       allocate(Methods(method_id)%halo_pointer(num_halos), stat=ierror)
00242       if (ierror /= 0) then
00243          ierror = prism_error_alloc
00244          call psmile_error (ierror, 'halo_pointer', (/num_halos/), 1, &
00245                             __FILE__, __LINE__ )
00246          return
00247       endif
00248 
00249       do i = 1, num_halos
00250 
00251          Methods(method_id)%halo_pointer(i)%halo_local = .false.
00252          Methods(method_id)%halo_pointer(i)%halo_shape = psmile_undef
00253          do j = 1, ndim_3d
00254             nullify(Methods(method_id)%halo_pointer(i)%halo_real(j)%vector)
00255             nullify(Methods(method_id)%halo_pointer(i)%halo_dble(j)%vector)
00256 #if defined ( PRISM_QUAD_TYPE )
00257             nullify(Methods(method_id)%halo_pointer(i)%halo_quad(j)%vector)
00258 #endif
00259          enddo ! j = 1, ndim_3d
00260 
00261       enddo ! i = 1, num_halos
00262 
00263    end subroutine
00264 
00265    ! ===========================================================================
00266 
00267    ! TODO: the last argument face, is not yet in use. It will be necessary when
00268    ! cyclic coordinate handling is implemented for halos, because theoretically
00269    ! a block could get two halos from another block. In this case just using
00270    ! grid ids would not be sufficient.
00271    ! The tag being generate is currently defined as follows:
00272    ! tag bits  1-16 : halotag
00273    ! tag bits 17-24 : local grid id of sender (assumes that the id is <= 255)
00274    ! tag bits 25-32 : local grid id of receiver (assums that the id is <= 255)
00275    ! this assumes that the tag has at least 32 bit
00276    integer function generate_send_tag (grid_id, remote_grid_id, face)
00277 
00278       integer, intent (in) :: grid_id
00279       integer, intent (in) :: remote_grid_id
00280       integer, intent (in) :: face
00281 
00282 #ifdef PRISM_ASSERTION
00283       if ((grid_id > 255) .or. (remote_grid_id > 255)) then
00284          call psmile_assert (__FILE__, __LINE__, &
00285                              "local grid id is bigger than 255 (halo tag definition needs to be adjusted)")
00286       endif
00287 #endif
00288       generate_send_tag = halotag + ishft(grid_id, 16) + ishft(remote_grid_id, 24)
00289 
00290    end function generate_send_tag
00291 
00292    ! ===========================================================================
00293 
00294    integer function generate_recv_tag (grid_id, remote_grid_id, face)
00295 
00296       integer, intent (in) :: grid_id
00297       integer, intent (in) :: remote_grid_id
00298       integer, intent (in) :: face
00299 
00300 #ifdef PRISM_ASSERTION
00301       if ((grid_id > 255) .or. (remote_grid_id > 255)) then
00302          call psmile_assert (__FILE__, __LINE__, &
00303                              "local grid id is bigger than 255 (halo tag definition needs to be adjusted)")
00304       endif
00305 #endif
00306       generate_recv_tag = halotag + ishft(grid_id, 24) + ishft(remote_grid_id, 16)
00307 
00308    end function generate_recv_tag
00309 
00310    ! ===========================================================================
00311 
00312    subroutine extract_data_2d_dble (tgt, tgt_shape, src, src_shape)
00313 
00314       integer, intent (in)           :: tgt_shape(2,ndim_2d)
00315 
00316       double precision, intent (out) :: tgt(tgt_shape(1,1):tgt_shape(2,1), 
00317                                             tgt_shape(1,2):tgt_shape(2,2))
00318 
00319       integer, intent (in)           :: src_shape(2,ndim_2d)
00320 
00321       double precision, intent (in)  :: src(src_shape(1,1):src_shape(2,1), 
00322                                             src_shape(1,2):src_shape(2,2))
00323 
00324       tgt = src(tgt_shape(1,1):tgt_shape(2,1), &
00325                 tgt_shape(1,2):tgt_shape(2,2))
00326 
00327    end subroutine extract_data_2d_dble
00328 
00329    ! ===========================================================================
00330 
00331    subroutine extract_data_2d_real (tgt, tgt_shape, src, src_shape)
00332 
00333       integer, intent (in) :: tgt_shape(2,ndim_2d)
00334 
00335       real, intent (out)   :: tgt(tgt_shape(1,1):tgt_shape(2,1), 
00336                                   tgt_shape(1,2):tgt_shape(2,2))
00337 
00338       integer, intent (in) :: src_shape(2,ndim_2d)
00339 
00340       real, intent (in)    :: src(src_shape(1,1):src_shape(2,1), 
00341                                   src_shape(1,2):src_shape(2,2))
00342 
00343       tgt = src(tgt_shape(1,1):tgt_shape(2,1), &
00344                 tgt_shape(1,2):tgt_shape(2,2))
00345 
00346    end subroutine extract_data_2d_real
00347 
00348    ! ===========================================================================
00349 
00350    subroutine extract_data_3d_dble (tgt, tgt_shape, src, src_shape)
00351 
00352       integer, intent (in)           :: tgt_shape(2,ndim_3d)
00353 
00354       double precision, intent (out) :: tgt(tgt_shape(1,1):tgt_shape(2,1), 
00355                                             tgt_shape(1,2):tgt_shape(2,2), 
00356                                             tgt_shape(1,3):tgt_shape(2,3))
00357 
00358       integer, intent (in)           :: src_shape(2,ndim_3d)
00359 
00360       double precision, intent (in)  :: src(src_shape(1,1):src_shape(2,1), 
00361                                            src_shape(1,2):src_shape(2,2), 
00362                                            src_shape(1,3):src_shape(2,3))
00363 
00364       tgt = src(tgt_shape(1,1):tgt_shape(2,1), &
00365                 tgt_shape(1,2):tgt_shape(2,2), &
00366                 tgt_shape(1,3):tgt_shape(2,3))
00367 
00368    end subroutine extract_data_3d_dble
00369 
00370    ! ===========================================================================
00371 
00372    subroutine extract_data_3d_real (tgt, tgt_shape, src, src_shape)
00373 
00374       integer, intent (in) :: tgt_shape(2,ndim_3d)
00375 
00376       real, intent (out)   :: tgt(tgt_shape(1,1):tgt_shape(2,1), 
00377                                   tgt_shape(1,2):tgt_shape(2,2), 
00378                                   tgt_shape(1,3):tgt_shape(2,3))
00379 
00380       integer, intent (in) :: src_shape(2,ndim_3d)
00381 
00382       real, intent (in)    :: src(src_shape(1,1):src_shape(2,1), 
00383                                   src_shape(1,2):src_shape(2,2), 
00384                                   src_shape(1,3):src_shape(2,3))
00385 
00386       tgt = src(tgt_shape(1,1):tgt_shape(2,1), &
00387                 tgt_shape(1,2):tgt_shape(2,2), &
00388                 tgt_shape(1,3):tgt_shape(2,3))
00389 
00390    end subroutine extract_data_3d_real
00391 
00392    ! ===========================================================================
00393 
00394    subroutine send_halo (method_id, halo, shift, ndim, ierror)
00395 
00396       integer, intent (in)                      :: method_id
00397       type (halo_info), intent (inout)          :: halo
00398       integer, intent (in)                      :: shift(ndim_3d)
00399       integer, intent (in)                      :: ndim
00400       integer, intent (inout)                   :: ierror
00401 
00402       integer                      :: curr_buffer_size
00403       integer                      :: local_send_range(2, ndim_3d)
00404       integer                      :: i
00405       integer                      :: tag
00406       type (coords_block), pointer :: coords_pointer
00407 
00408       tag = generate_send_tag(Methods(method_id)%grid_id, halo%remote_grid_id, 1)
00409 
00410       local_send_range(1,:) = halo%send_range(1,:) + shift
00411       local_send_range(2,:) = halo%send_range(2,:) + shift
00412 
00413 #ifdef DEBUG
00414       print 9995, trim(ch_id), " sending halo to rank ", halo%remote_rank, " with tag ", tag
00415       print 9996, trim(ch_id), "    local  range: ", local_send_range
00416       print 9996, trim(ch_id), "    global range: ", halo%send_range
00417       call psmile_flushstd
00418 #endif
00419 
00420       coords_pointer => Methods(method_id)%coords_pointer
00421 
00422       if (coords_pointer%coords_datatype == MPI_REAL) then
00423 
00424          ! check the buffer
00425          if (.not. associated (real_send_buf)) then
00426             curr_buffer_size = 0
00427          else
00428             curr_buffer_size = size (real_send_buf)
00429          endif
00430 
00431          ! increase size of buffer if necessary
00432          if (curr_buffer_size < halo%halo_size) then
00433             real_send_buf => psmile_realloc (real_send_buf, halo%halo_size)
00434          endif
00435 
00436          do i = 1, ndim
00437             ! copy halo into buffer
00438             if (ndim == ndim_3d) then
00439                call extract_data_3d_real (real_send_buf, local_send_range, &
00440                                           coords_pointer%coords_real(i)%vector, &
00441                                           coords_pointer%coords_shape)
00442             else
00443                call extract_data_2d_real (real_send_buf, local_send_range, &
00444                                           coords_pointer%coords_real(i)%vector, &
00445                                           coords_pointer%coords_shape)
00446             endif
00447 
00448             call psmile_bsend (real_send_buf, halo%halo_size, MPI_REAL, halo%remote_rank, &
00449                                tag, Comps(comp_id)%comm, ierror)
00450          enddo
00451 
00452       else ! (coords_pointer%coords_datatype == MPI_REAL)
00453 
00454          ! check the buffer
00455 
00456          if (.not. associated (dble_send_buf)) then
00457             curr_buffer_size = 0
00458          else
00459             curr_buffer_size = size (dble_send_buf)
00460          endif
00461 
00462          ! increase size of buffer if necessary
00463 
00464          if (curr_buffer_size < halo%halo_size) then
00465             dble_send_buf => psmile_realloc (dble_send_buf, halo%halo_size)
00466          endif
00467 
00468          do i = 1, ndim
00469 
00470             ! copy halo into buffer
00471 
00472             if (ndim == ndim_3d) then
00473                call extract_data_3d_dble (dble_send_buf, local_send_range, &
00474                                           coords_pointer%coords_dble(i)%vector, &
00475                                           coords_pointer%coords_shape)
00476             else
00477                call extract_data_2d_dble (dble_send_buf, local_send_range, &
00478                                           coords_pointer%coords_dble(i)%vector, &
00479                                           coords_pointer%coords_shape)
00480             endif
00481 
00482             call psmile_bsend (dble_send_buf, halo%halo_size, MPI_DOUBLE_PRECISION, halo%remote_rank, &
00483                                tag, Comps(comp_id)%comm, ierror)
00484          enddo
00485 
00486       endif ! (coords_pointer%coords_datatype == MPI_REAL)
00487 
00488 9995 format (1x, a, ': psmile_get_halo_points:', a, i4, a, i10)
00489 9996 format (1x, a, ': psmile_get_halo_points:', a, 6i4)
00490 
00491    end subroutine send_halo
00492 
00493    ! ===========================================================================
00494 
00495    subroutine recv_halo (method_id, halo, halo_id, shift, ndim, request, ierror)
00496 
00497       integer, intent (in) :: method_id
00498       type (halo_info), intent (inout) :: halo
00499       integer, intent (in) :: halo_id
00500       integer, intent (in) :: shift(ndim_3d)
00501       integer, intent (in) :: ndim
00502       integer, intent (out) :: request(ndim)
00503       integer, intent (inout) :: ierror
00504 
00505       integer                    :: tag
00506       integer                    :: i
00507       type (halo_block), pointer :: halo_pointer
00508 
00509       tag = generate_recv_tag (Methods(method_id)%grid_id, halo%remote_grid_id, 1)
00510 
00511       halo_pointer => Methods(method_id)%halo_pointer(halo_id)
00512 
00513       halo_pointer%halo_shape(1,:) = halo%global_range(1,:) + shift(:)
00514       halo_pointer%halo_shape(2,:) = halo%global_range(2,:) + shift(:)
00515 
00516 #ifdef DEBUG
00517       print 9995, trim(ch_id), " requesting halo from rank ", halo%remote_rank, " with tag ", tag
00518       print 9996, trim(ch_id), "    local  range: ", halo_pointer%halo_shape
00519       print 9996, trim(ch_id), "    global range: ", halo%global_range
00520       call psmile_flushstd
00521 #endif
00522 
00523       if (Methods(method_id)%coords_pointer%coords_datatype == MPI_REAL) then
00524 
00525          do i = 1, ndim
00526 
00527             allocate(halo_pointer%halo_real(i)%vector(halo%halo_size), &
00528                      stat=ierror)
00529             if (ierror /= 0) then
00530                ierror = PRISM_Error_Alloc
00531                call psmile_error ( ierror, 'halo_pointer', (/halo%halo_size/), 1, &
00532                      __FILE__, __LINE__ )
00533                return
00534             endif
00535 
00536             call MPI_Irecv (halo_pointer%halo_real(i)%vector, &
00537                             halo%halo_size, MPI_REAL, halo%remote_rank, tag, &
00538                             Comps(comp_id)%comm, request(i), ierror)
00539          enddo ! i = 1, ndim
00540 
00541       else
00542 
00543          do i = 1, ndim
00544 
00545             allocate(halo_pointer%halo_dble(i)%vector(halo%halo_size), &
00546                      stat=ierror)
00547             if (ierror /= 0) then
00548                ierror = PRISM_Error_Alloc
00549                call psmile_error ( ierror, 'halo_pointer', (/halo%halo_size/), 1, &
00550                      __FILE__, __LINE__ )
00551                return
00552             endif
00553 
00554             call MPI_Irecv (halo_pointer%halo_dble(i)%vector, &
00555                             halo%halo_size, MPI_DOUBLE_PRECISION, halo%remote_rank, tag, &
00556                             Comps(comp_id)%comm, request(i), ierror)
00557          enddo ! i = 1, ndim
00558 
00559       endif ! Methods(method_id)%coords_pointer%coords_datatype
00560 
00561 9995 format (1x, a, ': psmile_get_halo_points:', a, i4, a, i10)
00562 9996 format (1x, a, ': psmile_get_halo_points:', a, 6i4)
00563 
00564    end subroutine recv_halo
00565 
00566 end subroutine psmile_get_halo_points