Skip to content

Commit

Permalink
While determining masked blocks, take reentrancy and tripolar stitch …
Browse files Browse the repository at this point in the history 
…into account
  • Loading branch information
@alperaltuntas
alperaltuntas committed Nov 4, 2023
1 parent f2f4f04 commit f329b39
Show file tree
Hide file tree
Showing 2 changed files with 81 additions and 34 deletions.
13 changes: 12 additions & 1 deletion config_src/infra/FMS1/MOM_domain_infra.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ module MOM_domain_infra
public :: domain2D, domain1D, group_pass_type
! These interfaces are actually implemented or have explicit interfaces in this file.
public :: create_MOM_domain, clone_MOM_domain, get_domain_components, get_domain_extent
public :: deallocate_MOM_domain, get_global_shape, compute_block_extent
public :: deallocate_MOM_domain, get_global_shape, compute_block_extent, compute_extent
public :: pass_var, pass_vector, fill_symmetric_edges, rescale_comp_data
public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete
public :: create_group_pass, do_group_pass, start_group_pass, complete_group_pass
Expand Down Expand Up @@ -1945,6 +1945,17 @@ subroutine compute_block_extent(isg, ieg, ndivs, ibegin, iend)
call mpp_compute_block_extent(isg, ieg, ndivs, ibegin, iend)
end subroutine compute_block_extent

!> Get the array ranges in one dimension for the divisions of a global index space
subroutine compute_extent(isg, ieg, ndivs, ibegin, iend)
integer, intent(in) :: isg !< The starting index of the global index space
integer, intent(in) :: ieg !< The ending index of the global index space
integer, intent(in) :: ndivs !< The number of divisions
integer, dimension(:), intent(out) :: ibegin !< The starting index of each division
integer, dimension(:), intent(out) :: iend !< The ending index of each division

call mpp_compute_extent(isg, ieg, ndivs, ibegin, iend)
end subroutine compute_extent

!> Broadcast a 2-d domain from the root PE to the other PEs
subroutine broadcast_domain(domain)
type(domain2d), intent(inout) :: domain !< The domain2d type that will be shared across PEs.
Expand Down
102 changes: 69 additions & 33 deletions src/framework/MOM_domains.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -299,7 +299,8 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, &
! Auto-generate a mask file and determine the layout
call cpu_clock_begin(id_clock_auto_mask)
if (is_root_PE()) then
call gen_auto_mask_table(n_global, PEs_used, param_file, inputdir, auto_mask_table_fname, layout)
call gen_auto_mask_table(n_global, reentrant, tripolar_N, PEs_used, param_file, inputdir, &
auto_mask_table_fname, layout)
endif
call broadcast(layout, length=2)
call cpu_clock_end(id_clock_auto_mask)
Expand Down Expand Up @@ -434,38 +435,42 @@ subroutine MOM_define_layout(n_global, ndivs, layout)
end subroutine MOM_define_layout

!> Given a desired number of active npes, generate a layout and mask_table
subroutine gen_auto_mask_table(n_global, npes, param_file, inputdir, filename, layout)
subroutine gen_auto_mask_table(n_global, reentrant, tripolar_N, npes, param_file, inputdir, filename, layout)
integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions
integer, intent(in) :: npes !< The desired number of active PEs.
logical, dimension(2), intent(in) :: reentrant !< True if the x- and y- directions are periodic.
logical :: tripolar_N !< A flag indicating whether there is n. tripolar connectivity
integer, intent(in) :: npes !< The desired number of active PEs.
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
character(len=128), intent(in) :: inputdir !< INPUTDIR parameter>
character(len=:), allocatable, intent(in) :: filename !< Mask table file path (to be auto-generated.)
integer, dimension(2), intent(out) :: layout !< The generated layout of PEs (incl. masked blocks)
!local
real, dimension(n_global(1), n_global(2)) :: D ! Bathymetric depth (to be read in from TOPO_FILE)
integer, dimension(n_global(1), n_global(2)) :: mask ! Cell masks (based on D and MINIMUM_DEPTH)
real, dimension(n_global(1), n_global(2)) :: D ! Bathymetric depth (to be read in from TOPO_FILE)
integer, dimension(:,:), allocatable :: mask ! Cell masks (based on D and MINIMUM_DEPTH)
character(len=200) :: topo_filepath, topo_file ! Strings for file/path
character(len=200) :: topo_varname ! Variable name in file
character(len=200) :: topo_config
character(len=40) :: mdl = "gen_auto_mask_table" ! This subroutine's name.
integer :: i, j, p
real :: Dmask ! The depth for masking in the same units as D
real :: min_depth ! The minimum ocean depth in the same units as D
real :: mask_depth ! The depth shallower than which to mask a point as land.
real :: Dmask ! The depth for masking in the same units as D
real :: min_depth ! The minimum ocean depth in the same units as D
real :: mask_depth ! The depth shallower than which to mask a point as land.
real :: glob_ocn_frac ! ratio of ocean points to total number of points
real :: r_p ! aspect ratio for division count p.
real :: r_p ! aspect ratio for division count p.
integer :: nx, ny ! global domain sizes
integer, parameter :: ibuf=2, jbuf=2
real, parameter :: r_extreme = 10.0 ! aspect ratio limit (>1) for a layout to be considered.
integer :: num_masked_blocks
integer, allocatable :: mask_table(:,:)

! Read in params necessary for auto-masking
call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, do_not_log=.true., units="m", default=0.0)
call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, do_not_log=.true., units="m", default=0.0)
call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, do_not_log=.true., units="m", default=-9999.0)
call get_param(param_file, mdl, "TOPO_CONFIG", topo_config, do_not_log=.true., fail_if_missing=.true.)
call get_param(param_file, mdl, "TOPO_FILE", topo_file, do_not_log=.true., default="topog.nc")
call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, do_not_log=.true., default="depth")
topo_filepath = trim(inputdir)//trim(topo_file)

! Sanity checks
if (.not. is_root_pe()) then
call MOM_error(FATAL, 'gen_auto_mask_table should only be called by the root PE.')
Expand All @@ -477,36 +482,62 @@ subroutine gen_auto_mask_table(n_global, npes, param_file, inputdir, filename, l
call MOM_error(FATAL, " gen_auto_mask_table: Unable to open "//trim(topo_filepath))
endif

nx = n_global(1)
ny = n_global(2)

! Read in bathymetric depth.
D(:,:) = -9.0e30 ! Initializing to a very large negative depth (tall mountains) everywhere.
call read_field(topo_filepath, trim(topo_varname), D, start=(/1, 1/), nread=n_global, no_domain=.true.)

allocate( mask(nx+2*ibuf, ny+2*jbuf), source=0)

! Determine cell masks
Dmask = mask_depth
if (mask_depth == -9999.0) Dmask = min_depth
do i=1,n_global(1) ; do j=1,n_global(2)
do i=1,nx ; do j=1,ny
if (D(i,j) <= Dmask) then
mask(i,j) = 0
mask(i+ibuf,j+jbuf) = 0
else
mask(i,j) = 1
mask(i+ibuf,j+jbuf) = 1
endif
enddo ; enddo

glob_ocn_frac = real(sum(mask)) / size(mask)
! fill in buffer cells

if (reentrant(1)) then ! REENTRANT_X
mask(1:ibuf, :) = mask(nx+1:nx+ibuf, :)
mask(ibuf+nx+1:nx+2*ibuf, :) = mask(ibuf+1:2*ibuf, :)
endif

if (reentrant(2)) then ! REENTRANT_Y
mask(:, 1:jbuf) = mask(:, ny+1:ny+jbuf)
mask(:, jbuf+ny+1:ny+2*jbuf) = mask(:, jbuf+1:2*jbuf)
endif

if (tripolar_N) then ! TRIPOLAR_N
do i=1,nx+2*ibuf
do j=1,jbuf
mask(i, jbuf+ny+j) = mask(nx+2*ibuf+1-i, jbuf+ny+1-j)
enddo
enddo
!mask(:, jbuf+ny) = 1 ! TODO - REMOVE
endif

glob_ocn_frac = real(sum(mask(1+ibuf:nx+ibuf, 1+jbuf:ny+jbuf))) / (nx * ny)

! Iteratively check for all possible division counts starting from the upper bound of npes/glob_ocn_frac,
! which is over-optimistic for realistic domains, but may be satisfied with idealized domains.
! which is over-optimistic for realistic domains, but may be satisfied with idealized domains.
do p = ceiling(npes/glob_ocn_frac), npes, -1

! compute the layout for the current division count, p
call MOM_define_layout(n_global, p, layout)

! don't bother checking this p if the aspect ratio is extreme
r_p = (real(n_global(1))/layout(1)) / (real(n_global(2))/layout(2))
r_p = (real(nx)/layout(1)) / (real(ny)/layout(2))
if ( r_p * r_extreme < 1 .or. r_extreme < r_p ) cycle

! Get the number of masked_blocks for this particular division count
call determine_land_blocks(mask, n_global, layout(1), layout(2), num_masked_blocks)
call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks)

! If we can eliminate enough blocks to reach the target npes, adopt
! this p (and the associated layout) and terminate the iteration.
Expand All @@ -521,42 +552,47 @@ subroutine gen_auto_mask_table(n_global, npes, param_file, inputdir, filename, l
"of MOM6 PEs or set AUTO_MASKTABLE to False.")
endif

! Call determine_land_blocks once again, this time to retrieve the mask_table.
! Call determine_land_blocks once again, this time to retrieve and write out the mask_table.
allocate(mask_table(num_masked_blocks,2))
call determine_land_blocks(mask, n_global, layout(1), layout(2), num_masked_blocks, mask_table)
call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks, mask_table)
call write_auto_mask_file(mask_table, layout, npes, filename)
deallocate(mask_table)
deallocate(mask)

end subroutine gen_auto_mask_table

! Given number of domain divisions, compute the max number of land blocks that can be eliminated,
! Given a number of domain divisions, compute the max number of land blocks that can be eliminated,
! and return the resulting mask table if requested.
subroutine determine_land_blocks(mask, n_global, idiv, jdiv, num_masked_blocks, mask_table)
subroutine determine_land_blocks(mask, nx, ny, idiv, jdiv, ibuf, jbuf, num_masked_blocks, mask_table)
integer, dimension(:,:), intent(in) :: mask !< cell masks based on depth and MINIMUM_DEPTH
integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions
integer, intent(in) :: idiv !< number of divisions along x axis
integer, intent(in) :: jdiv !< number of divisions along y axis
integer, intent(in) :: nx !< Total number of gridpoints in x-dir (global)
integer, intent(in) :: ny !< Total number of gridpoints in y-dir (global)
integer, intent(in) :: idiv !< number of divisions along x-dir
integer, intent(in) :: jdiv !< number of divisions along y-dir
integer, intent(in) :: ibuf !< number of buffer cells in x-dir.
!! (not necessarily the same as NIHALO)
integer, intent(in) :: jbuf !< number of buffer cells in y-dir.
!! (not necessarily the same as NJHALO)
integer, intent(out) :: num_masked_blocks
integer, intent(out), optional :: mask_table(:,:)
! integer
integer, dimension(idiv) :: ibegin !< The starting index of each division along x axis
integer, dimension(idiv) :: iend !< The ending index of each division along x axis
integer, dimension(jdiv) :: jbegin !< The starting index of each division along y axis
integer, dimension(jdiv) :: jend !< The ending index of each division along y axis
integer, parameter :: ibuf=2, jbuf=2
integer :: i, j, ib, ie, jb,je

call compute_extent(1, n_global(1), idiv, ibegin, iend)
call compute_extent(1, n_global(2), jdiv, jbegin, jend)
call compute_extent(1, nx, idiv, ibegin, iend)
call compute_extent(1, ny, jdiv, jbegin, jend)

num_masked_blocks = 0

do i=1,idiv
ib = max(ibegin(i)-ibuf, 1)
ie = min(iend(i)+ibuf, n_global(1))
ib = ibegin(i)
ie = iend(i) + 2 * ibuf
do j=1,jdiv
jb = max(jbegin(j)-jbuf, 1)
je = min(jend(j)+jbuf, n_global(2))
jb = jbegin(j)
je = jend(j) + 2 * jbuf

if (any(mask(ib:ie,jb:je)==1)) cycle

Expand Down

0 comments on commit f329b39

Please sign in to comment.