CGpop mini-Application (2-sided MPI 1D data structure version) 0.1
simple_domain.F90
Go to the documentation of this file.
00001 !==============================================================================
00002 ! Copyright (C) 2010, University Corporation for Atmospheric Research,
00003 !                     Colorado State University,
00004 !                     Los Alamos National Security, LLC,
00005 !                     United States Department of Energy
00006 !
00007 ! All rights reserved.  See ../COPYING for copyright details
00008 !==============================================================================
00009 
00010 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
00011 !>
00012 !! Model domain and routines for initializing
00013 !! the domain.  This module also initializes the decompositions and
00014 !! distributions across processors/threads by calling relevent
00015 !! routines in the block, distribution modules.
00016 !<
00017 module simple_domain
00018     use kinds_mod, only: i4, r8, char_len, log_kind
00019     use simple_type, only: distrb
00020     use constants, only: blank_fmt, delim_fmt
00021     use communicate, only: MPI_COMM_OCN, my_task, master_task, get_num_procs
00022     use simple_blocks, only: nx_block,ny_block,nblocks_tot,nblocks_x, &
00023         nblocks_y, nghost, get_block_parameter, set_block_parameter, &
00024         all_blocks, all_blocks_ij, Neigh, GetEWBlockIndex, GetNSBlockIndex, &
00025         get_block_id, NumNeigh, ieast, iwest, isouth, inorth, iseast, iswest, &
00026         ineast, inwest
00027     use broadcast, only: broadcast_scalar
00028     use exit_mod, only: sigAbort, exit_pop
00029     use IOUnitsMod, only: stdout,nmlin,nml_filename
00030     use domain_size, only: nx_global,ny_global,max_blocks_tropic, &
00031         block_size_x,block_size_y
00032     use reductions, only: global_maxval
00033     use io_serial, only: read_AppMD, read_tiles
00034 
00035     implicit none
00036     private
00037     save
00038 
00039     ! !PUBLIC MEMBER FUNCTIONS
00040     public  :: init_domain_blocks ,&
00041         init_domain_distribution
00042 
00043     public :: read_solverioT
00044     public :: read_local_tiles 
00045 
00046     !>
00047     !! Reads blocks or tiles that are specific for each task. 
00048     !<
00049     interface read_local_tiles
00050         module procedure read_local_tiles_int
00051         module procedure read_local_tiles_dbl
00052     end interface
00053 
00054     ! !PUBLIC DATA MEMBERS:
00055     integer (i4), public :: 
00056         nblocks_tropic     ! actual number of blocks in tropic distribution
00057 
00058     integer (i4), dimension(:), pointer, public :: 
00059         blocks_tropic      ! block ids for local blocks in barotropic dist
00060 
00061     !------------------------------------------------------------
00062     type (distrb), public ::  !  block distribution info
00063         distrb_tropic      ! block distribution for barotropic part
00064 
00065     !-----------------------------------------------------------------------
00066     !
00067     !   module private variables - for the most part these appear as
00068     !   module variables to facilitate sharing info between init_domain1
00069     !   and init_domain2.
00070     !
00071     !-----------------------------------------------------------------------
00072     character (char_len) ::      
00073         ew_boundary_type,        ! type of domain bndy in each logical
00074         ns_boundary_type           !    direction (ew is i, ns is j)
00075 
00076     integer (i4) ::     ! decomposition info
00077         nprocs_tropic    ! num of processors in barotropic dist
00078 
00079     !***********************************************************************
00080     contains
00081 
00082     !>
00083     !!  This routine reads in domain information and calls the routine
00084     !!  to set up the block decomposition.
00085     !!
00086     !! @param reorder An integer array that reorders block from the cannonical 
00087     !!                 order in the input files to a space-filling curve order 
00088     !!            which is used by the miniapp.
00089     !<
00090     subroutine init_domain_blocks(reorder)
00091         integer (i4), pointer :: reorder(:)
00092 
00093         !------------------------------------------------------------------
00094         !
00095         !  local variables
00096         !
00097         !------------------------------------------------------------------
00098 
00099         integer (i4) :: 
00100             nml_error          ! namelist read error flag
00101 
00102         integer (i4) :: i,j,n
00103 
00104         integer (i4) :: iblock_src,jblock_src
00105         integer (i4) :: iblock_east, jblock_east, 
00106             iblock_west, jblock_west, 
00107             iblock_south,jblock_south, 
00108             iblock_north,jblock_north, 
00109             iblock_swest,jblock_swest, 
00110             iblock_nwest,jblock_nwest, 
00111             iblock_seast,jblock_seast, 
00112             iblock_neast,jblock_neast
00113 
00114         character(len=80) :: tilefile
00115 
00116         !------------------------------------------------------------------
00117         !
00118         !  input namelists
00119         !
00120         !------------------------------------------------------------------
00121         namelist /input_nml/ tilefile
00122 
00123         !------------------------------------------------------------------
00124         !
00125         !  read domain information from namelist input
00126         !
00127         !------------------------------------------------------------------
00128         nprocs_tropic = get_num_procs()
00129         ew_boundary_type = 'cyclic'
00130         ns_boundary_type = 'closed'
00131 
00132         !------------------------------------------------------------------
00133         !
00134         !  perform some basic checks on domain
00135         !
00136         !------------------------------------------------------------------
00137         if (nx_global < 1 .or. ny_global < 1) then
00138             !***
00139             !*** domain size zero or negative
00140             !***
00141             call exit_POP(sigAbort,'Invalid domain: size < 1') ! no domain
00142         else if (nghost < 2) then
00143             !***
00144             !*** must have at least 2 layers of ghost cells
00145             !***
00146             call exit_POP(sigAbort,'Not enough ghost cells allocated')
00147         endif
00148 
00149         !------------------------------------------------------------------
00150         !
00151         !  compute block decomposition and details
00152         !
00153         !------------------------------------------------------------------
00154 
00155         if (my_task == master_task) then
00156         open (nmlin, file=nml_filename, status='old',iostat=nml_error)
00157         if (nml_error /= 0) then
00158         nml_error = -1
00159         else
00160         nml_error =  1
00161         endif
00162         do while (nml_error > 0)
00163         read(nmlin, nml=input_nml,iostat=nml_error)
00164         end do
00165         if (nml_error == 0) close(nmlin)
00166         endif
00167         call broadcast_scalar(tilefile,       master_task)
00168 
00169         call read_AppMD(tilefile,all_blocks,reorder)
00170         nblocks_x   = (nx_global-1)/block_size_x + 1
00171         nblocks_y   = (ny_global-1)/block_size_y + 1
00172         nblocks_tot = nblocks_x*nblocks_y
00173 
00174         allocate(all_blocks_ij(nblocks_x,nblocks_y))
00175         n=0
00176         do j=1,nblocks_y
00177         do i=1,nblocks_x
00178         n = n + 1
00179         all_blocks_ij(i,j) = n
00180         enddo
00181         enddo
00182 
00183         allocate(Neigh(NumNeigh,nblocks_tot))
00184 
00185 
00186         do n=1,nblocks_tot
00187         call get_block_parameter(n,iblock=iblock_src,jblock=jblock_src)
00188 
00189         !*** compute cartesian i,j block indices for each neighbor
00190         !*** use zero if off the end of closed boundary
00191         !*** use jnorth=nblocks_y and inorth < 0 for tripole boundary
00192         !***   to make sure top boundary communicated to all top
00193         !***   boundary blocks
00194 
00195         call GetEWBlockIndex(ew_boundary_type,iblock_src,jblock_src, &
00196         iblock_east,jblock_east,iblock_west,jblock_west)
00197 
00198         call GetNSBlockIndex(ns_boundary_type,iblock_src,jblock_src, &
00199         iblock_north,jblock_north,iblock_south,jblock_south)
00200 
00201         call GetEWBlockIndex(ew_boundary_type,iblock_south,jblock_south, &
00202         iblock_seast,jblock_seast,iblock_swest,jblock_swest)
00203 
00204         call GetEWBlockIndex(ew_boundary_type,iblock_north,jblock_north, &
00205         iblock_neast,jblock_neast,iblock_nwest,jblock_nwest)
00206 
00207         ! save Neighbor information
00208         Neigh(ieast,n)  = get_block_id(iblock_east,jblock_east)
00209         Neigh(iwest,n)  = get_block_id(iblock_west,jblock_west)
00210         Neigh(inorth,n) = get_block_id(iblock_north,jblock_north)
00211         Neigh(isouth,n) = get_block_id(iblock_south,jblock_south)
00212         Neigh(iseast,n) = get_block_id(iblock_seast,jblock_seast)
00213         Neigh(ineast,n) = get_block_id(iblock_neast,jblock_neast)
00214         Neigh(iswest,n) = get_block_id(iblock_swest,jblock_seast)
00215         Neigh(inwest,n) = get_block_id(iblock_nwest,jblock_nwest)
00216         enddo
00217 
00218         !------------------------------------------------------------------
00219         !
00220         !  Now we need grid info before proceeding further
00221         !  Print some domain information
00222         !
00223         !------------------------------------------------------------------
00224 
00225         if (my_task == master_task) then
00226         write(stdout,delim_fmt)
00227         write(stdout,blank_fmt)
00228         write(stdout,'(a18)') 'Domain Information'
00229         write(stdout,blank_fmt)
00230         write(stdout,delim_fmt)
00231         write(stdout,'(a26,i6)') '  Horizontal domain: nx = ',nx_global
00232         write(stdout,'(a26,i6)') '                     ny = ',ny_global
00233         write(stdout,'(a26,i6)') '  Block size:  nx_block = ',nx_block
00234         write(stdout,'(a26,i6)') '               ny_block = ',ny_block
00235         write(stdout,'(a29,i6)') '  Processors for barotropic: ', &
00236         nprocs_tropic
00237         write(stdout,'(a25,i2)') '  Number of ghost cells: ', nghost
00238         endif
00239     end subroutine init_domain_blocks
00240 
00241     !>
00242     !! This routine calls appropriate setup routines to distribute blocks
00243     !! across processors and defines arrays with block ids for any local
00244     !! blocks. Information about ghost cell update routines is also
00245     !! initialized here through calls to the appropriate boundary routines.
00246     !<
00247     subroutine init_domain_distribution(reorder)
00248         ! !INPUT PARAMETERS:
00249         integer (i4),pointer,  intent(inout) :: reorder(:)
00250         integer(i4), save :: timer_distrib
00251 
00252         !------------------------------------------------------------------
00253         !
00254         !  local variables
00255         !
00256         !------------------------------------------------------------------
00257         character (char_len) :: outstring
00258 
00259         integer (i4), parameter :: 
00260         max_work_unit=10   ! quantize the work into values from 1,max
00261 
00262         integer (i4) :: 
00263             i,j,k,n              ,! dummy loop indices
00264             count1, count2       ,! dummy counters
00265             work_unit            ,! size of quantized work unit
00266             nblocks_tmp          ,! temporary value of nblocks
00267             nblocks_tmp_tropic   ,! num blocks on proc for tropic
00268             nblocks_max_tropic     ! max blocks on proc for tropic
00269 
00270         integer (i4), dimension(:), allocatable :: 
00271             nocn               ,! number of ocean points per block
00272             work_per_block       ! number of work units per block
00273 
00274         integer (i4), dimension(:), allocatable :: 
00275             noWork_per_block     ! =1 for all land blocks
00276 
00277         integer (i4), dimension(:), pointer :: iglob,jglob
00278         integer (i4) :: jb,je, ib, ie
00279 
00280         type (distrb) :: distrb_calc
00281         integer(i4), pointer :: Lreorder(:)
00282 
00283         !------------------------------------------------------------------
00284         !  estimate the amount of work per processor using the topography
00285         !------------------------------------------------------------------
00286         allocate(nocn(nblocks_tot))
00287         nocn = 0
00288         nocn(:) = all_blocks(:)%npoints
00289 
00290         !------------------------------------------------------------------
00291         !
00292         !  determine the distribution of blocks across processors
00293         !
00294         !------------------------------------------------------------------
00295         distrb_tropic = calc_distribution(nprocs_tropic,nocn,reorder)
00296         deallocate(nocn)
00297 
00298         !------------------------------------------------------------------
00299         !
00300         !  allocate and determine block id for any local blocks in each
00301         !  distribution.
00302         !
00303         !------------------------------------------------------------------
00304         call create_local_block_ids(blocks_tropic, distrb_tropic)
00305 
00306         if (my_task < distrb_tropic%nprocs .and. &
00307             associated(blocks_tropic)) then
00308             nblocks_tropic = size(blocks_tropic)
00309         else
00310             nblocks_tropic = 0
00311         endif
00312 
00313         nblocks_max_tropic = global_maxval(nblocks_tropic)
00314 
00315         if (nblocks_max_tropic > max_blocks_tropic) then
00316             write(outstring,*) 'tropic blocks exceed max: increase max from ', &
00317             max_blocks_tropic, ' to ',nblocks_max_tropic
00318             call exit_POP(sigAbort,trim(outstring))
00319         else if (nblocks_max_tropic < max_blocks_tropic) then
00320             write(outstring,*) 'tropic blocks too large: decrease max to',&
00321             nblocks_max_tropic
00322             if (my_task == master_task) write(stdout,*) trim(outstring)
00323         endif
00324     end subroutine init_domain_distribution
00325 
00326     !>
00327     !!  This subroutine partitions the linear spacing-filling curve into 
00328     !!  a number of equal length segments.
00329     !!
00330     !! @param nprocs    The number of tasks over which to partition the problem.
00331     !! @param nocn      An integer array which indicates the number of active
00332     !!                  ocean points within a block.
00333     !! @param reorder   An integer array which reorders the blocks from
00334     !!                  cannonical to space-filling curve order.
00335     !<
00336     function  calc_distribution(nprocs,nocn,reorder) result(dist)
00337         integer (i4), intent(in) :: nprocs
00338         integer (i4), intent(in) :: nocn(:)
00339         integer (i4), intent(in) :: reorder(:)
00340         type (distrb) :: dist
00341 
00342         integer (i4) :: nb,nblocks,nActive,nblocksL,extra,s1
00343         integer (i4) :: n,il,nblocksT,pid
00344         integer (i4) :: ii,tmp1
00345         integer (i4), allocatable :: proc_tmp(:)
00346 
00347         nActive = COUNT(reorder(:) > 0) 
00348         nb = SIZE(nocn)
00349 
00350         dist%nprocs = nprocs
00351         dist%communicator = MPI_COMM_OCN 
00352         allocate(dist%proc(nb))
00353         allocate(dist%local_block(nb))
00354 
00355         nblocks=nActive  ! This is the total number of active blocks
00356         nblocksL = nblocks/nprocs
00357 
00358         ! every cpu gets nblocksL blocks, but the first 'extra' get nblocksL+1
00359         extra = mod(nblocks,nprocs)
00360         s1 = extra*(nblocksL+1)
00361 
00362         dist%proc(:)=0
00363         dist%local_block(:)=0
00364 
00365         !---------------------------------------
00366         ! use the cryptic SFC decomposition code 
00367         !---------------------------------------
00368         do n=1,nblocks_tot
00369         ii = reorder(n)
00370         if(ii>0) then 
00371         if(ii<=s1) then 
00372         ii=ii-1
00373         tmp1 = ii/(nblocksL+1)
00374         dist%proc(n) = tmp1+1
00375         else
00376         ii=ii-s1-1
00377         tmp1 = ii/nblocksL
00378         dist%proc(n) = extra + tmp1 + 1
00379         endif
00380         endif
00381         enddo
00382 
00383         allocate(proc_tmp(nprocs))
00384         proc_tmp = 0
00385         do n=1,nblocks_tot
00386         pid = dist%proc(n)
00387         if(pid>0) then
00388         proc_tmp(pid) = proc_tmp(pid) + 1
00389         dist%local_block(n) = proc_tmp(pid)
00390         endif
00391         enddo
00392         deallocate(proc_tmp)
00393     end function calc_distribution
00394 
00395     !>
00396     !! This routine determines which blocks in an input distribution are
00397     !! located on the local processor and creates an array of block ids
00398     !! for all local blocks.
00399     !<
00400     subroutine create_local_block_ids(block_ids, distribution)
00401         ! !INPUT PARAMETERS:
00402         type (distrb), intent(in) :: 
00403         distribution           ! input distribution for which local
00404 
00405         ! !OUTPUT PARAMETERS:
00406         integer (i4), dimension(:), pointer :: 
00407         block_ids              ! array of block ids for every block
00408                                ! that resides on the local processor
00409 
00410         !-------------------------------------------------------------------
00411         !
00412         !  local variables
00413         !
00414         !--------------------------------------------------------------------
00415         integer (i4) :: 
00416         n, bid, bcount        ! dummy counters
00417 
00418         !--------------------------------------------------------------------
00419         !
00420         !  first determine number of local blocks to allocate array
00421         !
00422         !--------------------------------------------------------------------
00423         bcount = 0
00424         do n=1,size(distribution%proc)
00425             if (distribution%proc(n) == my_task+1) bcount = bcount + 1
00426         end do
00427 
00428         if (bcount > 0) allocate(block_ids(bcount))
00429 
00430         !-------------------------------------------------------------------
00431         !
00432         !  now fill array with proper block ids
00433         !
00434         !-------------------------------------------------------------------
00435         if (bcount > 0) then
00436             do n=1,size(distribution%proc)
00437                 if (distribution%proc(n) == my_task+1) then
00438                     block_ids(distribution%local_block(n)) = n
00439                 endif
00440             end do
00441         endif
00442     end subroutine create_local_block_ids
00443 
00444     !>
00445     !! This subroutine reads in the blocks or tiles assigned to each task.
00446     !! Note that for simplicity this subroutine calls read_tiles which 
00447     !! reads in all blocks.  Blocks not owned by the task are discarded.
00448     !!
00449     !! @param fname   The name of the netCDF from which to read the blocks.
00450     !! @param vname   The name of the variable to read from the netCDF file.
00451     !! @param var     The variable into which data will be read.
00452     !<
00453     subroutine read_local_tiles_dbl(fname,vname,var)
00454         character(len=80), intent(in) :: fname
00455         character(len=*), intent(in) :: vname
00456         real(r8), dimension(:,:,:) :: var
00457         real(r8), allocatable :: tmp(:,:,:)
00458         integer(i4) :: n, gbid
00459 
00460         allocate(tmp(nx_block,ny_block,nblocks_tot))
00461         call read_tiles(fname,vname,tmp)
00462 
00463         ! just pull out tiles that
00464         do n=1,nblocks_tropic
00465             gbid = blocks_tropic(n)
00466             var(:,:,n) = tmp(:,:,gbid)
00467         enddo
00468         deallocate(tmp)
00469     end subroutine read_local_tiles_dbl
00470 
00471     !>
00472     !! This subroutine reads in the blocks or tiles assigned to each task.
00473     !! Note that for simplicity this subroutine calls read_tiles which 
00474     !! reads in all blocks.  Blocks not owned by the task are discarded.
00475     !!
00476     !! @param fname   The name of the netCDF from which to read the blocks.
00477     !! @param vname   The name of the variable to read from the netCDF file.
00478     !! @param var     The variable into which data will be read.
00479     !<
00480     subroutine read_local_tiles_int(fname,vname,var)
00481         character(len=80), intent(in) :: fname
00482         character(len=*), intent(in) :: vname
00483         integer(i4), dimension(:,:,:) :: var
00484         real(r8), allocatable :: tmp(:,:,:)
00485         integer(i4) :: n, gbid
00486 
00487         allocate(tmp(nx_block,ny_block,nblocks_tot))
00488         call read_tiles(fname,vname,tmp)
00489 
00490         ! just pull out tiles that
00491         do n=1,nblocks_tropic
00492             gbid = blocks_tropic(n)
00493             var(:,:,n) = tmp(:,:,gbid)
00494         enddo
00495         deallocate(tmp)
00496     end subroutine read_local_tiles_int
00497 
00498     !>
00499     !! This subroutine reads in all solver state necessary to initialize 
00500     !! the miniapp.
00501     !!
00502     !! @param fname    The filename of the netCDF tilefile which contains the
00503     !!                 solver state.
00504     !! @param nstep    The timestep in the simulation that the inputfile
00505     !!                 represents.
00506     !! @param A0       The diagonal coefficient of the stencil which forms the
00507     !!                 sparse matrix vector multiply.
00508     !! @param AN       The northern neighbor coefficient of the stencil.
00509     !! @param AE       The eastern neighbor coefficient of the stencil.
00510     !! @param ANE      The northeastern neighbor coefficient of the stencil.
00511     !! @param RHS      The right and side of the linear system. This
00512     !!                 corresponds to 'b' from the equation 'Ax=b'.
00513     !! @param PRESSI   The initial guess for the surface pressure.  This
00514     !!                 corresponds to 'x'.
00515     !! @param PRESSF   The final estimate for the surface pressure.
00516     !! @param RCALCT   The land mask.  Values greater than 0 indicate ocean
00517     !!                 points.
00518     !! @param TAREASQ  The area for each grid point.
00519     !! @param KMT      This integer array which indicates the number of
00520     !!                 vertical levels.
00521     !<
00522     subroutine read_solverioT(fname,nstep,A0,AN,AE,ANE, &
00523         RHS,PRESSI,PRESSF,RCALCT, TAREASQ,KMT)
00524 
00525         character(len=80), intent(in) :: fname
00526         integer(i4), intent(inout) :: nstep
00527         real(r8), intent(inout), 
00528             dimension(nx_block,ny_block,max_blocks_tropic) :: A0, AN, AE, ANE
00529         real(r8), intent(inout), 
00530             dimension(nx_block,ny_block,max_blocks_tropic) :: PRESSI,RHS,PRESSF
00531         real(r8), intent(inout), 
00532             dimension(nx_block,ny_block,max_blocks_tropic) :: RCALCT,TAREASQ
00533         integer(i4), intent(inout), optional, 
00534             dimension(nx_block,ny_block,max_blocks_tropic) :: KMT
00535 
00536         call read_local_tiles(fname,'A0',A0)
00537         call read_local_tiles(fname,'AN',AN)
00538         call read_local_tiles(fname,'AE',AE)
00539         call read_local_tiles(fname,'ANE',ANE)
00540         call read_local_tiles(fname,'PRESS_I',PRESSI)
00541         call read_local_tiles(fname,'RHS',RHS)
00542         call read_local_tiles(fname,'PRESS_F',PRESSF)
00543         call read_local_tiles(fname,'RCALCT',RCALCT)
00544         call read_local_tiles(fname,'TAREASQ',TAREASQ)
00545 
00546         if(present(KMT)) then
00547             call read_local_tiles(fname,'KMT',KMT)
00548         endif
00549     end subroutine read_solverioT
00550 end module simple_domain
00551 
00552 !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 All Classes Namespaces Files Functions Variables