Skip to content

Updates to li_flood_fill #144

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

Updates to li_flood_fill #144

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9d7d456
Updates to li_flood_fill
matthewhoffman Apr 24, 2025
2a4086e
Sum masked cells over nCellsSolve instead of nCells
matthewhoffman Apr 29, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
136 changes: 55 additions & 81 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1233,7 +1233,7 @@ subroutine remove_small_islands(domain, err)
end where

call mpas_log_write("***Cleaning up stranded cells after removing small islands***")
call li_flood_fill(seedMask, growMask, domain)
call li_flood_fill(domain)
do iCell = 1, nCells
if (li_mask_is_floating_ice(cellMask(iCell)) .and. seedMask(iCell) == 0) then
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell)
Expand Down Expand Up @@ -3272,7 +3272,7 @@ subroutine find_open_ocean(domain, openOceanMask, err)
call mpas_dmpar_field_halo_exch(domain, 'seedMask')
call mpas_timer_stop("halo updates")

call li_flood_fill(seedMask, growMask, domain)
call li_flood_fill(domain)
openOceanMask = seedMask

call mpas_deallocate_scratch_field(seedMaskField, single_block_in=.true.)
Expand Down Expand Up @@ -4001,7 +4001,7 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d
growMask = 1
end where

call li_flood_fill(seedMask, growMask, domain)
call li_flood_fill(domain)

! Remove ice from flood-filled mask, and any neighboring non-dynamic ice
do iCell = 1, nCells
Expand Down Expand Up @@ -4307,7 +4307,7 @@ subroutine li_remove_icebergs(domain)
where ( (seedMask == 0) .and. li_mask_is_floating_ice(cellMask(:)) .and. li_mask_is_dynamic_ice(cellMask(:)) )
growMask = 1
endwhere
call li_flood_fill(seedMask, growMask, domain)
call li_flood_fill(domain)

! Add floating non-dynamic fringe, but exclude dynamic ice isolated by
! non-dynamic ice, which can cause velocity solver to fail to converge.
Expand All @@ -4319,7 +4319,7 @@ subroutine li_remove_icebergs(domain)
elsewhere
growMask = 0
endwhere
call li_flood_fill(seedMask, growMask, domain)
call li_flood_fill(domain)

! Now remove any ice that was not flood-filled - these are icebergs
block => domain % blocklist
Expand Down Expand Up @@ -4372,25 +4372,29 @@ end subroutine li_remove_icebergs
!> The calling routine defines the seedMask to specify the initial masked region, and the growMask
!> to specify the potential region to be filled by the flood-fill routine. The flood-fill routine
!> then adds the contiguous areas of the growMask to the seedMask, which is passed back to the
!> calling routine for application. Any integer fields in Registry can be used for seedMask and
!> growMask, but there are existing variables with those names that can be used.
!> calling routine for application.
!> The calling routine must use the fields seedMask and growMask in the scratch pool in Registry.
!> This is because the halo updates in this routine are hard-coded to use those field names.
!> Because those two fields are scratch variables, the calling routine needs to be responsible
!> for allocating (and populating) those fields prior to calling this routine, and then deallocating
!> them afterwards.
!> \author Trevor Hillebrand and Matthew Hoffman
!> \date March 2021
!> \details
!-----------------------------------------------------------------------
subroutine li_flood_fill(seedMask, growMask, domain)
subroutine li_flood_fill(domain)
!-----------------------------------------------------------------
! input/output variables
!-----------------------------------------------------------------
type (domain_type), intent(inout) :: domain !< Input/Output: domain object
integer, dimension(:), intent(inout) :: seedMask
integer, dimension(:), intent(in) :: growMask

! Local variables
type (block_type), pointer :: block
type (mpas_pool_type), pointer :: meshPool
type (mpas_pool_type), pointer :: geometryPool
type (mpas_pool_type), pointer :: velocityPool
type (mpas_pool_type), pointer :: scratchPool

integer, dimension(:), pointer :: seedMask
integer, dimension(:), pointer :: growMask
integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell
integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell
integer, dimension(:), allocatable :: seedMaskOld
Expand All @@ -4404,41 +4408,21 @@ subroutine li_flood_fill(seedMask, growMask, domain)
err = 0

block => domain % blocklist
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(scratchPool, 'seedMask', seedMask)
call mpas_pool_get_array(scratchPool, 'growMask', growMask)

allocate(seedMaskOld(nCells+1))
seedMaskOld(:) = 0

!call mpas_log_write("Flood-fill: allocated.")

! First mark grounded ice to initialize flood fill mask
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)

! make sure masks are up to date. May not be necessary, but safer to
! do anyway.
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
err = ior(err, err_tmp)
newMaskCountLocal = sum(seedMask(1:nCellsSolve))

!call mpas_log_write("Flood-fill: updated masks.")
newMaskCountLocal = sum(seedMask)
!call mpas_log_write("Initialized $i cells to local mask", intArgs=(/newMaskCountLocal/))

block => block % next
end do

!call mpas_log_write("Flood-fill initialization complete.")

! Outer loop over processors (should also have a loop over blocks)
! Inner loop over cells on that processor
! Outer loop over processors; Inner loop over cells on that processor

! Initialize global mask count
call mpas_dmpar_sum_int(domain % dminfo, newMaskCountLocal, newMaskCountGlobal)
Expand All @@ -4460,46 +4444,36 @@ subroutine li_flood_fill(seedMask, growMask, domain)

! Now update (advance) mask locally

block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)

! initialize local loop
localLoopCount = 0
newMaskCountLocal = 1 ! need to make sure we enter the loop
do while (newMaskCountLocal > 0)
localLoopCount = localLoopCount + 1
!call mpas_log_write(" Starting local cell loop $i", intArgs=(/localLoopCount/))

! initialize
newMaskCountLocal = 0
seedMaskOld(:) = seedMask(:)

do iCell = 1, nCellsSolve ! this gives owned cells only
if ( growMask(iCell)>0 ) then
! If it has a marked neighbor, then add it to the mask
do n = 1, nEdgesOnCell(iCell)
jCell = cellsOnCell(n, iCell)
if ( (seedMaskOld(jCell) == 1) .and. (seedMask(iCell) .ne. 1) ) then
seedMask(iCell) = 1
newMaskCountLocal = newMaskCountLocal + 1
exit ! skip the rest of this do-loop - no need to check additional neighbors
endif
enddo
endif ! if not already marked
enddo ! loop over cells

! Accumulate cells added locally until we do the next global
! reduce
newMaskCountLocalAccum = newMaskCountLocalAccum + newMaskCountLocal
!call mpas_log_write(" Added $i new cells to local mask", intArgs=(/newMaskCountLocal/))
enddo ! local mask loop

block => block % next
end do
! initialize local loop
localLoopCount = 0
newMaskCountLocal = 1 ! need to make sure we enter the loop
do while (newMaskCountLocal > 0)
localLoopCount = localLoopCount + 1
!call mpas_log_write(" Starting local cell loop $i", intArgs=(/localLoopCount/))

! initialize
newMaskCountLocal = 0
seedMaskOld(:) = seedMask(:)

do iCell = 1, nCellsSolve ! this gives owned cells only
if (growMask(iCell) > 0) then
! If it has a marked neighbor, then add it to the mask
do n = 1, nEdgesOnCell(iCell)
jCell = cellsOnCell(n, iCell)
if ( (seedMaskOld(jCell) == 1) .and. (seedMask(iCell) .ne. 1) ) then
seedMask(iCell) = 1
newMaskCountLocal = newMaskCountLocal + 1
exit ! skip the rest of this do-loop - no need to check additional neighbors
endif
enddo
endif ! if not already marked
enddo ! loop over cells

! Accumulate cells added locally until we do the next global
! reduce
newMaskCountLocalAccum = newMaskCountLocalAccum + newMaskCountLocal
!call mpas_log_write(" Added $i new cells to local mask", intArgs=(/newMaskCountLocal/))
enddo ! local mask loop

! update count of cells added to mask globally
call mpas_dmpar_sum_int(domain % dminfo, newMaskCountLocalAccum, newMaskCountGlobal)
Expand Down
2 changes: 1 addition & 1 deletion components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -531,7 +531,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
growMask(iCell) = 1
end if
end do
call li_flood_fill(seedMask, growMask, domain)
call li_flood_fill(domain)

floatingBasalMassBal(:) = floatingBasalMassBal(:) * seedMask(:)

Expand Down
2 changes: 1 addition & 1 deletion components/mpas-albany-landice/src/mode_forward/mpas_li_ocean_extrap.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ subroutine li_ocean_extrap_solve(domain, err)
endif
enddo
! now create mask of all marine locations connected to open ocean - to be used below to screen out lakes
call li_flood_fill(connectedMarineMask, growMask, domain)
call li_flood_fill(domain)

! make it a 3D mask based on the topography (loop through nISMIP6OceanLayers)
call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers)
Expand Down