DOUG 0.2

Distribution_elem.F90

Go to the documentation of this file.
00001 module Distribution_elem_mod
00002   use Vect_mod
00003   use Mesh_plot_mod
00004   use ElemMtxs_mods
00005   use SpMtx_util
00006   use Distribution_base_mod
00007 
00008   implicit none
00009 
00010 #include<doug_config.h>
00011 
00012 ! "on-the-fly" real/complex picking
00013 #ifdef D_COMPLEX
00014 #define float complex
00015 #else
00016 #define float real
00017 #endif
00018 
00019   private
00020   public :: parallelAssembleFromElemInput, Distribution_elem_addoverlap
00021 
00022 contains
00023 
00024   !----------------------------------------------------------------
00027   subroutine parallelAssembleFromElemInput(Msh, A, &
00028                b, nparts, part_opts, A_interf)
00029     implicit none
00030 
00031     type(Mesh),     intent(in out) :: Msh !< Mesh
00032     type(SpMtx),    intent(out) :: A !< System matrix
00033     float(kind=rk), dimension(:), pointer :: b !< local RHS
00034     ! Partitioning
00035     integer, intent(in) :: nparts !< number of parts to partition a mesh
00036     integer, dimension(6), intent(in) :: part_opts !< partition options (see METIS manual)
00037     type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface
00038 
00039     if (ismaster()) then ! MASTER
00040        write(stream,*)
00041        write(stream,*) 'master thread'
00042        if (D_MSGLVL > 1) &
00043             call MasterCtrlData_print()
00044 
00045     else ! SLAVES
00046        write(stream,'(a,i4,a)') 'slave [',myrank,'] thread'
00047        if (D_MSGLVL > 1) &
00048             call SharedCtrlData_print()
00049     end if
00050 
00051     ! =======================
00052     ! Mesh and its Dual Graph
00053     ! =======================
00054     !
00055     ! Create Mesh object
00056 
00057     Msh = Mesh_New()
00058 
00059     if (ismaster()) then
00060        ! Initialise Mesh object
00061        call Mesh_initFromFile(Msh, trim(mctls%info_file))
00062     endif
00063 
00064     ! Get from master Mesh's parameters: nell, ngf, mfrelt, nsd, nnode
00065     call Mesh_paramsMPIBCAST(Msh)
00066     if (D_MSGLVL > 1) &
00067          call Mesh_printInfo(Msh)
00068 
00069     ! Master reads in from files: feedom lists, coordinates, freedom map
00070     if (ismaster()) then
00071        call Mesh_readFromFile(Msh, &
00072             fnFreelists = trim(mctls%freedom_lists_file), &
00073             fnCoords    = trim(mctls%coords_file),        &
00074             fnFreemap   = trim(mctls%freemap_file))
00075     end if
00076 
00077     ! Master non-blockingly sends mesh data: nfrelt, mhead
00078     call Mesh_dataMPIISENDRECV(Msh, &
00079          nfrelt   = .true., &
00080          mhead    = .true.)
00081     if (D_MSGLVL > 4) &
00082          call Mesh_printElemFree(Msh)
00083 
00084     ! For multi-variable problems which have more than one block
00085     if (sctls%number_of_blocks > 1) then
00086        if (ismaster()) then
00087           call Mesh_readFromFile(Msh, &
00088                fnFreemask = trim(mctls%freedom_mask_file))
00089        end if
00090        call Mesh_dataMPIISENDRECV(Msh, freemask = .true.)
00091     end if
00092 
00093 
00094     ! Build dual graph (Graph object is a data field in Mesh class)
00095     ! let all procs do this - compare whith broadcasting the one bult on master
00096     call Mesh_buildGraphDual(Msh)
00097 
00098     ! Partition mesh's dual graph
00099     if (ismaster()) then
00100        if (sctls%plotting == D_PLOT_YES) then
00101           call Mesh_pl2D_mesh(Msh)
00102        end if
00103 
00104        ! Partition graph: D_PART_PMETIS, D_PART_KMETIS, D_PART_VKMETIS
00105        call Mesh_partitionDual(Msh, nparts, D_PART_VKMETIS, part_opts)
00106        if (sctls%plotting == D_PLOT_YES) then
00107           call Mesh_pl2D_partitions(Msh)
00108        end if
00109     endif
00110 
00111     ! Distribute elements to partitons map among slaves
00112     call Mesh_dataMPIISENDRECV(Msh, eptnmap=.true.)
00113 
00114     ! Calculate number of elements in partitions
00115     call Mesh_calcElemsInParts(Msh)
00116 
00117     ! Build global to local, local to global maps and
00118     ! inner/interface mask for freedoms
00119     ! (also finds local number of freedoms 'Mesh%nlf';
00120     !  this hidden appearance of 'Mesh_findNLF()' helps
00121     !  to speed up calculations a bit)
00122     call Mesh_buldMapsNghbrsMasksNLF(Msh)
00123 
00124     ! ===============================
00125     ! Assemble system matrix and RHS
00126     ! ===============================
00127     A = SpMtx_New()
00128     allocate(b(Msh%nlf))
00129     b = 0.0_rk
00130 
00131     if (ismaster()) then
00132        if (numprocs>1.and.present(A_interf)) then
00133           A_interf = SpMtx_New()
00134           call ElemMtxs_readAndDistribute(Msh, trim(mctls%elemmat_rhs_file), A, b, A_interf)
00135        else
00136           call ElemMtxs_readAndDistribute(Msh, trim(mctls%elemmat_rhs_file), A, b)
00137        end if
00138     else
00139        if (numprocs>1.and.present(A_interf)) then
00140           A_interf = SpMtx_New()
00141           call ElemMtxs_recvAndAssemble(Msh, A, b, A_interf)
00142        else
00143           call ElemMtxs_recvAndAssemble(Msh, A, b)
00144        end if
00145     end if
00146 
00147     if (sctls%verbose>9) then
00148       write(stream,'(/a)') 'System matrix:'
00149       call SpMtx_printInfo(A)
00150       if (A%nrows <= 25) then
00151          call SpMtx_printMat(A)
00152       else if ((A%nrows > 25).and.(A%nrows <= 100)) then
00153          call SpMtx_printRaw(A)
00154       end if
00155     endif
00156 
00157     ! ==================
00158     ! Finish assemble local RHS
00159     ! ==================
00160     if (A%nrows <= 25) & ! if (D_MSGLVL > 2) &
00161          call Vect_Print(b, 'RHS assembled (local) ')
00162 
00163     ! initialise auxiliary data for manipulating vectors
00164     call Vect_setIntfEnd(sum(Msh%inner_interf_fmask))
00165     call Vect_buildDotMask(Msh)
00166     if (D_MSGLVL > 4) &
00167          call Vect_Print(dot_intf_fmask,'dot_intf_fmask ')
00168     call Vect_buildDotMap()
00169     if (D_MSGLVL > 4) &
00170          call Vect_Print(dot_intf_fmap,'dot_intf_fmap ')
00171 
00172     ! Free mesh graph, not needed anymore
00173     call Mesh_destroyGraph(Msh)
00174 
00175     ! Update inner node count
00176     Msh%ninner=Msh%nlf
00177 
00178   end subroutine parallelAssembleFromElemInput
00179 
00180   subroutine Distribution_elem_addoverlap(D,x)
00181     type(Distribution),intent(in) :: D
00182     float(kind=rk),dimension(:),intent(in out)   :: x ! Vector
00183     float(kind=rk),dimension(:),pointer          :: x_tmp ! TMP Vector
00184     integer :: i, n, p, count
00185     ! MPI
00186     integer, dimension(:), pointer :: in_reqs
00187     integer                        :: ierr, out_req, status(MPI_STATUS_SIZE)
00188     integer, parameter             :: D_TAG_FREE_INTERFFREE = 777
00189 
00190     ! initialise receives
00191     allocate(in_reqs(D%mesh%nnghbrs))
00192     allocate(x_tmp(size(x))) !TODO: remove this -- should not be needed
00193     x_tmp=x
00194     do p = 1,D%mesh%nparts
00195        if (D%mesh%nfreesend_map(p) /= 0) then
00196           i = D%cache%pid2indx(p)
00197           n = D%mesh%nfreesend_map(p)
00198           call MPI_IRECV(D%cache%inbufs(i)%arr, n, MPI_fkind, &
00199                p-1, D_TAG_FREE_INTERFFREE, MPI_COMM_WORLD, in_reqs(i), ierr)
00200        end if
00201     end do
00202     ! Need a barrier?
00203     call MPI_Barrier(MPI_COMM_WORLD,ierr)
00204     ! nonblockingly send
00205     do p = 1,D%mesh%nparts
00206        if (D%mesh%nfreesend_map(p) /= 0) then
00207           i = D%cache%pid2indx(p)
00208           n = D%mesh%nfreesend_map(p)
00209           D%cache%outbufs(i)%arr(1:n) = x_tmp(D%cache%fexchindx(1:n,i))
00210           call MPI_ISEND(D%cache%outbufs(i)%arr, n, MPI_fkind, &
00211                p-1, D_TAG_FREE_INTERFFREE, MPI_COMM_WORLD, out_req, ierr)
00212        end if
00213     end do
00214     ! Need a barrier?
00215     call MPI_Barrier(MPI_COMM_WORLD,ierr)
00216     !
00217     ! ...some work could be done here...
00218     !
00219     ! wait for neighbours' interface freedoms
00220     do while (.true.)
00221        call MPI_WAITANY(D%mesh%nnghbrs, in_reqs, i, status, ierr)
00222        if (i /= MPI_UNDEFINED) then
00223           count = D%mesh%nfreesend_map(D%mesh%nghbrs(i)+1)
00224           x(D%cache%fexchindx(1:count,i)) = &
00225                x(D%cache%fexchindx(1:count,i)) + D%cache%inbufs(i)%arr
00226        else
00227           exit
00228        end if
00229     end do
00230     deallocate(x_tmp)
00231     deallocate(in_reqs)
00232   end subroutine Distribution_elem_addoverlap
00233   
00234 end module Distribution_elem_mod