psmile_get_field_dble.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_field_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_get_field_dble (field_id, data_array, len, nbr_fields, &
00012                                         ierror)
00013 !
00014 ! !USES:
00015 !
00016       use PRISM_constants
00017       use PSMILe, dummy_interface => PSMILe_Get_field_dble
00018       use PSMILe_SMIOC, only : sga_smioc_comp,transient_in
00019 
00020       implicit none
00021 !
00022 ! !INPUT PARAMETERS:
00023 !
00024       Integer, Intent (In)             :: field_id
00025 
00026 !     Description of the field
00027 
00028       Integer, Intent (In)             :: len
00029 
00030 !     1st dimension of data_array
00031 
00032       Integer, Intent (In)             :: nbr_fields
00033 
00034 !     Number of fields to be sent
00035 
00036 ! !INPUT/OUTPUT PARAMETERS:
00037 !
00038       Double Precision, Intent (InOut) :: data_array (len, nbr_fields)
00039 !
00040 !     Field data
00041 
00042 ! !OUTPUT PARAMETERS:
00043 !
00044       Integer, Intent (Out)            :: ierror
00045 
00046 !     Returns the error code of PSMILe_Get_field_dble;
00047 !             ierror = 0 : No error
00048 !             ierror > 0 : Severe error
00049 !
00050 ! !LOCAL VARIABLES
00051 !
00052 !     ... field pointer
00053 !
00054       Type (Gridfunction), Pointer    :: field
00055       Type (Taskin_type), Pointer     :: fieldin
00056       Type (transient_in), Pointer    :: transi_in
00057 
00058       Logical                         :: conservation
00059 
00060       Integer                         :: trans_in_id
00061 !
00062 !     ... Method Id and pointer
00063 !
00064       Integer                         :: task_id
00065       Integer                         :: method_id
00066       Type (Method), Pointer          :: mp
00067       Integer                         :: var_shape (2, max_dim)
00068 !
00069 !     ... Grid information
00070 !
00071       Integer                         :: grid_id, grid_type
00072 !
00073 !     ... for direct communication
00074 !
00075       Integer                         :: status (MPI_STATUS_SIZE)
00076 !
00077       Integer                         :: index, n, tag, tag0
00078 !
00079       Integer, parameter              :: nd_msgdata = 3
00080       Integer                         :: msgdata (nd_msgdata)
00081 !
00082 !     ... for communication with the coupler process
00083 !
00084 #ifdef SEND_TO_COUPLER
00085 rename PRISM_INIT, ... in prismtrs.inc; name conflict
00086       Integer                         :: msgtrans (PRISM_Header_length)
00087 #endif
00088       type pbuff
00089          Double Precision, Pointer    :: p (:)
00090       end type pbuff
00091       type(pbuff), allocatable        :: pbuffer (:)
00092       double precision, allocatable   :: buffer (:)
00093 
00094 !     for global conservation
00095       Double Complex                  :: old_global_sum(nbr_fields)
00096       Double Complex                  :: new_global_sum(nbr_fields)
00097       Double Complex                  :: temp_global_sum(nbr_fields)
00098       Double Complex                  :: correction_factor(nbr_fields)
00099       Double Complex, Pointer         :: sum_buffer(:)
00100       Integer                         :: i, j
00101 !
00102 !     ... for error handling
00103 !
00104       Integer, parameter              :: nerrp = 3
00105       Integer                         :: ierrp (nerrp)
00106 !
00107 ! !DESCRIPTION:
00108 !
00109 ! Subroutine "PSMILe_Get_field_dble" receives the data (returned in
00110 ! data_array) from the sources of field "field".
00111 !
00112 !
00113 ! !REVISION HISTORY:
00114 !
00115 !   Date      Programmer   Description
00116 ! ----------  ----------   -----------
00117 ! 03.07.21    H. Ritzdorf  created
00118 !
00119 !EOP
00120 !----------------------------------------------------------------------
00121 !
00122 !  $Id: psmile_get_field_dble.F90 2755 2010-11-19 15:19:52Z hanke $
00123 !  $Author: hanke $
00124 !
00125    Character(len=len_cvs_string), save :: mycvs = 
00126        '$Id: psmile_get_field_dble.F90 2755 2010-11-19 15:19:52Z hanke $'
00127 !
00128 !----------------------------------------------------------------------
00129 !
00130 !  Initialization
00131 !
00132 #ifdef VERBOSE
00133       print 9990, trim(ch_id), field_id, trim(Fields (field_id)%local_name)
00134 
00135       call psmile_flushstd
00136 #endif /* VERBOSE */
00137 !
00138       ierror = 0
00139 !
00140 !===> Get field pointer
00141 !
00142       field => Fields (field_id)
00143       fieldin => Fields (field_id)%Taskin
00144       transi_in => sga_smioc_comp(field%comp_id)%sga_smioc_transi(field%smioc_loc)%sg_transi_in
00145 
00146 !     get conservation status from origin fields
00147 #ifdef PRISM_ASSERTION
00148       if (transi_in%ig_nb_in_orig .ne. 1) then
00149          call psmile_assert (__FILE__, __LINE__, "Multiple origins")
00150       endif
00151 #endif
00152       conservation = transi_in%sga_in_orig(1) %ig_conserv .ne. PSMILe_undef
00153 
00154 #ifdef PRISM_ASSERTION
00155       if ( field%status == PSMILe_Status_free .or. &
00156            field%status == PSMILe_Status_undefined ) then
00157          call psmile_assert (__FILE__, __LINE__, "Incorrect field")
00158       endif
00159 #endif
00160 !
00161 !===> Get method id and pointer to the method
00162 !
00163       method_id = field%method_id
00164       grid_id   = Methods(method_id)%grid_id
00165       grid_type = Grids(grid_id)%grid_type
00166 !
00167       mp => Methods(method_id)
00168 !
00169 #ifdef PRISM_ASSERTION
00170       if ( mp%status == PSMILe_Status_free .or. &
00171            mp%status == PSMILe_Status_undefined ) then
00172          call psmile_assert (__FILE__, __LINE__, "Incorrect method")
00173       endif
00174 #endif
00175 !
00176 !===> Receive the data which were directly sent by application processes
00177 !     (i)  Receive info concerning the data sent
00178 !     (ii) Receive data
00179 !
00180 !rr
00181 !rr Currently this only works for one input channel. We have to loop over task_ids
00182 !rr  in future.
00183 !rr
00184       task_id = 1
00185       tag0 = datatag + fieldin%In_Channel(task_id)%global_transi_id
00186 !
00187 #ifdef DEBUG
00188       print 9970, trim(ch_id), fieldin%n_recv_direct, fieldin%n_recv_coupler
00189 9970 format (1x, a, ': psmile_get_field_dble: #recvs: direct, coupler', 2i4)
00190 #endif
00191 !
00192       if ( grid_type == PRISM_Gaussreduced_regvrt ) then
00193          var_shape (:,1) = field%var_shape(:,1)
00194          var_shape (:,2) = 1
00195          var_shape (:,3) = field%var_shape(:,3)
00196          var_shape (:,4:max_dim) = field%var_shape(:,3:max_dim-1)
00197       else
00198          var_shape = field%var_shape
00199       endif
00200 !
00201       do n = 1, fieldin%n_recv_direct
00202 !
00203          index = fieldin%recv_direct(n)%recv_info_index
00204 !
00205 !        ... Receive info on data
00206 !
00207 ! msgdata (1) = Method id in destination process
00208 ! msgdata (2) = Datatype sent
00209 !               Hier schon auf kleineren Typ kopieren,
00210 !               wenn Empfaenger was kleineres braucht.
00211 !               Dazu brauchen wir aber die XML Infos.
00212 ! msgdata (3) = Number of fields sent
00213 !
00214 ! Dieses Receive ist unschon; MPI_ANY_SOURCE waere schoener
00215 ! aber dafuer muessen wird erst den Tag und die Liste eindeutig machen.
00216 !
00217          tag = tag0 + method_id*1000
00218 !      print *, psmile_rank, 'recv from ', mp%recv_infos_direct(index)%source, tag, method_id
00219 ! Hier gibt es einen Deadlock, wenn die Daten nicht gesendet werden
00220 ! Beispielsweise, wenn eine Funktion nicht gesendet wird.
00221 ! Kann man da eine vernuenftige Fehlermeldung generieren und
00222 ! einen Test einbauen ?
00223 !
00224 !        print *, 'Waiting for header: ', nd_msgdata, ' data from ', &
00225 !                       mp%recv_infos_direct(index)%source, ' with tag ', tag 
00226          call MPI_Recv (msgdata, nd_msgdata, MPI_INTEGER, &
00227                         mp%recv_infos_direct(index)%source, tag, comm_psmile, &
00228                         status, ierror)
00229 !
00230          if (ierror /= MPI_SUCCESS) then
00231             ierrp (1) = ierror
00232             ierrp (2) = status (MPI_SOURCE)
00233             ierrp (3) = status (MPI_TAG)
00234 
00235             ierror = PRISM_Error_Recv
00236 
00237             call psmile_error (ierror, 'MPI_Recv', &
00238                                ierrp, 3, __FILE__, __LINE__ )
00239             return
00240          endif
00241 !
00242          if (method_id /= msgdata (1)) then
00243             ierrp (1) = field_id
00244             ierrp (2) = method_id
00245             ierrp (3) = msgdata (1)
00246 
00247             ierror = PRISM_Error_Internal
00248 
00249             call psmile_error (ierror, 'Wrong message received', &
00250                                ierrp, 3, __FILE__, __LINE__ )
00251             return
00252             
00253          endif
00254 
00255          if (nbr_fields /= msgdata (3)) then
00256             ierrp (1) = field_id
00257             ierrp (2) = nbr_fields
00258             ierrp (3) = msgdata (3)
00259 
00260             ierror = PRISM_Error_Internal
00261 
00262             call psmile_error (ierror, 'Inconsistend number of fields', &
00263                                ierrp, 3, __FILE__, __LINE__ )
00264             return
00265             
00266          endif
00267 
00268          if (msgdata (2) /= PRISM_DOUBLE_PRECISION) then
00269             ierrp (1) = msgdata (2)
00270 
00271             ierror = PRISM_Error_Internal
00272 
00273             call psmile_error (ierror, 'Sorry: Different datatypes currently not supported', &
00274                                ierrp, 1, __FILE__, __LINE__ )
00275             return
00276          endif
00277 
00278 !        ... Receive the data
00279 
00280          call psmile_get_irr_field_dble (data_array, var_shape, nbr_fields, &
00281                                    mp%recv_infos_direct(index)%dstijk,  &
00282                                    mp%recv_infos_direct(index)%npoints, &
00283                                    mp%recv_infos_direct(index)%dstars,  &
00284                                    mp%recv_infos_direct(index)%nar,     &
00285                                    mp%recv_infos_direct(index)%nloc,    &
00286                                    status(MPI_SOURCE),                  &
00287                                    tag, comm_psmile, ierror)
00288          if (ierror > 0) return
00289 
00290       end do
00291 !
00292 !===> Receive the data which is sent by the coupler process
00293 !
00294 !     (i)  Receive request to the transformer
00295 !     (ii) Receive corresponding data to the transformer
00296 !
00297 #ifdef VERBOSE
00298       print 9971, trim(ch_id), conservation
00299 #endif
00300       if (.not. conservation) then
00301          do n = 1, fieldin%n_recv_coupler
00302 !
00303             index       = fieldin%recv_coupler(n)%recv_info_index
00304             trans_in_id = fieldin%recv_coupler(n)%trans_in_id
00305 !
00306 !===> ... Receive data
00307 !     PSMILe_Trs_get_dble aufbrechen !!!
00308 !
00309             Allocate (buffer(mp%recv_infos_coupler(index)%nloc*nbr_fields), stat = ierror)
00310             if ( ierror > 0 ) then
00311                ierrp (1) = ierror
00312                ierrp (2) = mp%recv_infos_coupler(index)%nloc * nbr_fields
00313 
00314                ierror = PRISM_Error_Alloc
00315                call psmile_error ( ierror, 'buffer', &
00316                                  ierrp, 2, __FILE__, __LINE__ )
00317                return
00318             endif
00319 
00320             call psmile_trs_get_dble (trans_in_id, &
00321                                     mp%recv_infos_coupler(index)%epio_id,         &
00322                                     mp%recv_infos_coupler(index)%trs_rank,        &
00323                                     mp%recv_infos_coupler(index)%nloc*nbr_fields, &
00324                                     buffer, nbr_fields, ierror)
00325             if (ierror > 0) return
00326    !
00327    !        Copy data from compact list into destination array
00328    !
00329             call psmile_put_compact_list_3d_dble (buffer, &
00330                                     mp%recv_infos_coupler(index)%dstijk,  &
00331                                     mp%recv_infos_coupler(index)%nloc,    &
00332                                     data_array, var_shape, nbr_fields,    &
00333                                     ierror)
00334             if (ierror > 0) return
00335    !
00336             Deallocate (buffer)
00337 
00338          end do
00339       else !.not. conservation
00340 
00341 !        allocate buffer pointers for all data packages
00342          allocate (pbuffer(fieldin%n_recv_coupler), stat = ierror)
00343          if ( ierror > 0 ) then
00344             ierrp (1) = ierror
00345             ierrp (2) = mp%recv_infos_coupler(index)%nloc * nbr_fields
00346 
00347             ierror = PRISM_Error_Alloc
00348             call psmile_error ( ierror, 'pbuffer', &
00349                               ierrp, 2, __FILE__, __LINE__ )
00350             return
00351          endif
00352 
00353          new_global_sum = 0
00354 
00355          do n = 1, fieldin%n_recv_coupler
00356 !
00357             index       = fieldin%recv_coupler(n)%recv_info_index
00358             trans_in_id = fieldin%recv_coupler(n)%trans_in_id
00359 !
00360 !===> ... Receive data
00361 !     PSMILe_Trs_get_dble aufbrechen !!!
00362 !
00363             Allocate (pbuffer(n)%p(mp%recv_infos_coupler(index)%nloc*nbr_fields), stat = ierror)
00364             if ( ierror > 0 ) then
00365                ierrp (1) = ierror
00366                ierrp (2) = mp%recv_infos_coupler(index)%nloc * nbr_fields
00367 
00368                ierror = PRISM_Error_Alloc
00369                call psmile_error ( ierror, 'buffer', &
00370                                  ierrp, 2, __FILE__, __LINE__ )
00371                return
00372             endif
00373 
00374             call psmile_trs_get_dble (trans_in_id, &
00375                                     mp%recv_infos_coupler(index)%epio_id,         &
00376                                     mp%recv_infos_coupler(index)%trs_rank,        &
00377                                     mp%recv_infos_coupler(index)%nloc*nbr_fields, &
00378                                     pbuffer(n)%p, nbr_fields, ierror)
00379             if (ierror > 0) return
00380 
00381 !           allocate buffer for computing the sum of the field
00382             Allocate (sum_buffer(mp%recv_infos_coupler(index)%nloc*nbr_fields), stat = ierror)
00383             if ( ierror > 0 ) then
00384                ierrp (1) = ierror
00385                ierrp (2) = mp%recv_infos_coupler(index)%nloc * nbr_fields
00386 
00387                ierror = PRISM_Error_Alloc
00388                call psmile_error ( ierror, 'sum_buffer', &
00389                                  ierrp, 2, __FILE__, __LINE__ )
00390                return
00391             endif
00392             ! initilise buffer of local data
00393             Where (pbuffer(n)%p == PSMILe_undef)
00394                sum_buffer = 0
00395             elsewhere
00396                sum_buffer = pbuffer(n)%p
00397             end where
00398             ! compute the global sum of the received field
00399             call psmile_global_sum_compute_dble(sum_buffer,&
00400                   mp%recv_infos_coupler(index)%nloc, nbr_fields,&
00401                   Comps(field%comp_id)%act_comm, temp_global_sum,&
00402                   ierror)
00403 
00404             do i = 1, nbr_fields
00405                call psmile_ddadd(temp_global_sum(i), new_global_sum(i))
00406             end do
00407             ! get the global sum of the original field
00408             call psmile_global_sum_recv_dble (old_global_sum, nbr_fields,&
00409                            mp%recv_infos_coupler(index)%trs_rank, ierror)
00410          end do
00411 
00412          ! compute the correction factor
00413          do i = 1, nbr_fields
00414             if (new_global_sum(i) .eq. 0) then
00415                correction_factor(i) = 0
00416             else
00417                correction_factor(i) = old_global_sum(i) / new_global_sum(i)
00418             endif
00419          end do
00420 #ifdef VERBOSE
00421          print 9972, trim(ch_id), 'new_global_sum=', new_global_sum
00422          print 9972, trim(ch_id), 'old_global_sum=', old_global_sum
00423          print 9972, trim(ch_id), 'correction_factor=', correction_factor
00424 #endif
00425          do n = 1, fieldin%n_recv_coupler
00426             index = fieldin%recv_coupler(n)%recv_info_index
00427             ! apply the correction factor
00428             do i = 1, nbr_fields
00429                do j = (i-1)*mp%recv_infos_coupler(index)%nloc+1, i*mp%recv_infos_coupler(index)%nloc
00430                   if (pbuffer(n)%p(j) .ne. PSMILe_undef) then
00431                      pbuffer(n)%p(j) = pbuffer(n)%p(j) * correction_factor(i)
00432                   endif
00433                end do
00434             end do
00435    !
00436    !        Copy data from compact list into destination array
00437    !
00438             call psmile_put_compact_list_3d_dble (pbuffer(n)%p, &
00439                                     mp%recv_infos_coupler(index)%dstijk,  &
00440                                     mp%recv_infos_coupler(index)%nloc,    &
00441                                     data_array, var_shape, nbr_fields,    &
00442                                     ierror)
00443             if (ierror > 0) return
00444          end do
00445 
00446 !        free allocated memory
00447          do n = 1, fieldin%n_recv_coupler
00448             Deallocate (pbuffer(n)%p)
00449          end do
00450          Deallocate (pbuffer)
00451       endif !.not. conservation
00452 !
00453 !===> All done
00454 !
00455 #ifdef VERBOSE
00456       print 9980, trim(ch_id), ierror
00457 
00458       call psmile_flushstd
00459 #endif /* VERBOSE */
00460 !
00461 !  Formats:
00462 !
00463 
00464 #ifdef VERBOSE
00465 
00466 9990 format (1x, a, ': psmile_get_field_dble: field_id', i5, '; name ', a)
00467 9980 format (1x, a, ': psmile_get_field_dble: eof ierror =', i3)
00468 9971 format (1x, a, ': psmile_get_field_dble: conservation status: ', L1)
00469 9972 format (1x, a, ': psmile_get_field_dble: ', a , 2(1X,e14.6))
00470 
00471 #endif /* VERBOSE */
00472 
00473       end subroutine PSMILe_get_field_dble