Skip to content

Do not update MPAS masks during advection #115

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

Merged
matthewhoffman merged 3 commits into from
Aug 26, 2024
Merged

Do not update MPAS masks during advection #115

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a975404
Do not update cell, edge, or vertex masks during RK loop
trhille Apr 12, 2024
f135bdf
Remove cellMask from vertical_remap
trhille Apr 12, 2024
f6a6653
Update masks before RK loop, but not at start of each advection stage.
trhille Aug 19, 2024
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
8 changes: 8 additions & 0 deletions components/mpas-albany-landice/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1433,6 +1433,14 @@ is the value of that variable from the *previous* time level!
description="temporary copy of cellMask"
persistence="scratch"
/>
<var name="edgeMaskTemporary" type="integer" dimensions="nEdges" units="none"
description="temporary copy of edgeMask"
persistence="scratch"
/>
<var name="vertexMaskTemporary" type="integer" dimensions="nVertices" units="none"
description="temporary copy of vertexMask"
persistence="scratch"
/>
</var_struct>

<!-- ================ -->
Expand Down
45 changes: 32 additions & 13 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,8 @@ subroutine li_advection_thickness_tracers(&

type (field1DInteger), pointer :: &
cellMaskTemporaryField, & ! scratch field containing old values of cellMask
edgeMaskTemporaryField, &
vertexMaskTemporaryField, &
thermalCellMaskField

! Allocatable arrays need for flux-corrected transport advection
Expand All @@ -201,6 +203,7 @@ subroutine li_advection_thickness_tracers(&
integer, dimension(:), pointer :: &
cellMask, & ! integer bitmask for cells
edgeMask, & ! integer bitmask for edges
vertexMask, &
thermalCellMask

integer, dimension(:,:), pointer :: cellsOnEdge
Expand Down Expand Up @@ -290,6 +293,7 @@ subroutine li_advection_thickness_tracers(&
call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1)
call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening)
call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)

Expand Down Expand Up @@ -356,6 +360,12 @@ subroutine li_advection_thickness_tracers(&
call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField)
call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.)

call mpas_pool_get_field(geometryPool, 'edgeMaskTemporary', edgeMaskTemporaryField)
call mpas_allocate_scratch_field(edgeMaskTemporaryField, .true.)

call mpas_pool_get_field(geometryPool, 'vertexMaskTemporary', vertexMaskTemporaryField)
call mpas_allocate_scratch_field(vertexMaskTemporaryField, .true.)

call mpas_pool_get_field(scratchPool, 'iceCellMask', thermalCellMaskField)
call mpas_allocate_scratch_field(thermalCellMaskField, .true.)
thermalCellMask => thermalCellMaskField % array
Expand All @@ -366,12 +376,14 @@ subroutine li_advection_thickness_tracers(&
! given the old thickness, compute the thickness in each layer
call li_calculate_layerThickness(meshPool, thickness, layerThickness)

! update masks
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
err = ior(err, err_tmp)

! save old copycellMask for determining cells changing from grounded to floating and vice versa
! Save copies of masks because we need to preserve mask
! states prior to advection for accurate time integration.
! A mask update is necessary to calculate grounding line flux,
Comment thread
matthewhoffman marked this conversation as resolved.
! after which we will reset the masks to their previous states.
cellMaskTemporaryField % array(:) = cellMask(:)
edgeMaskTemporaryField % array(:) = edgeMask(:)
vertexMaskTemporaryField % array(:) = vertexMask(:)

layerThicknessEdgeFlux(:,:) = 0.0_RKIND

!-----------------------------------------------------------------
Expand Down Expand Up @@ -540,10 +552,10 @@ subroutine li_advection_thickness_tracers(&
dynamicThickening = (sum(layerThickness, 1) - thickness) / dt * scyr ! units of m/yr


! Update the thickness and cellMask before applying the mass balance.
! The update is needed because the SMB and BMB depend on whether ice is present.
! Update the thickness before applying the mass balance, but
! do not update masks because mass balance acts on geometry
! before advection took place.
thickness = sum(layerThickness, 1)
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)



Expand Down Expand Up @@ -627,6 +639,9 @@ subroutine li_advection_thickness_tracers(&
enddo
endif

! We need an updated set of masks to calculate fluxAcrossGroundingLine,
Copy link
Copy Markdown

@mperego mperego Apr 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why the flux at the grounding line is calculated at each stage? Also, if I understand correctly at this point the thickness is not consistent with the velocity. Wouldn't be better to move the calculation of the GL flux outside of this function?

! but we will reset this to the previous state below for accuracy of the
! time integration scheme.
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
err = ior(err, err_tmp)

Expand Down Expand Up @@ -667,12 +682,17 @@ subroutine li_advection_thickness_tracers(&
endif
enddo ! edges

! Reset masks to state before advection and mass balance for
! accuracy of time integration scheme.
cellMask(:) = cellMaskTemporaryField % array(:)
edgeMask(:) = edgeMaskTemporaryField % array(:)
vertexMask(:) = vertexMaskTemporaryField % array(:)

! Remap tracers to the standard vertical sigma coordinate
! Note: If tracers are not being advected, then this subroutine simply restores the
! layer thickness to sigma coordinate values.

call vertical_remap(thickness, cellMask, meshPool, layerThickness, advectedTracers, err_tmp)
call vertical_remap(thickness, meshPool, layerThickness, advectedTracers, err_tmp)
err = ior(err, err_tmp)

if (config_print_thickness_advection_info) then
Expand Down Expand Up @@ -727,6 +747,8 @@ subroutine li_advection_thickness_tracers(&
call mpas_deallocate_scratch_field(basalTracersField, .true.)
call mpas_deallocate_scratch_field(surfaceTracersField, .true.)
call mpas_deallocate_scratch_field(cellMaskTemporaryField, .true.)
call mpas_deallocate_scratch_field(edgeMaskTemporaryField, .true.)
call mpas_deallocate_scratch_field(vertexMaskTemporaryField, .true.)
call mpas_deallocate_scratch_field(thermalCellMaskField, .true.)

! === error check
Expand Down Expand Up @@ -1939,7 +1961,7 @@ end subroutine li_layer_normal_velocity
!> OpenMP over either blocks or cells.
!
!-----------------------------------------------------------------------
subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers, err)
subroutine vertical_remap(thickness, meshPool, layerThickness, tracers, err)

!-----------------------------------------------------------------
!
Expand All @@ -1953,9 +1975,6 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers
real(kind=RKIND), dimension(:), intent(in) :: &
thickness !< Input: ice thickness

integer, dimension(:), intent(in) :: &
cellMask !< Input: mask for cells (needed for determining presence/absence of ice)

!-----------------------------------------------------------------
!
! input/output variables
Expand Down
5 changes: 2 additions & 3 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 @@ -3828,9 +3828,8 @@ subroutine li_finalize_damage_after_advection(domain, err)
call mpas_pool_get_array(geometryPool, 'damageNye', damageNye)
call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)

! 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)
! Note: In order to preserve accuracy of time integration,
! we do not update masks before finalizing damage.

if (config_preserve_damage) then
do iCell = 1, nCells
Expand Down
15 changes: 9 additions & 6 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -250,8 +250,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
err = ior(err, err_tmp)
call mpas_timer_stop("face melting for grounded ice")

! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? ***

! === Basal melting for floating ice ===========
call mpas_timer_start("basal melting for floating ice")
call li_basal_melt_floating_ice(domain, err_tmp)
Expand Down Expand Up @@ -364,7 +362,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
// ' is not supported with config_rk_order = $i', &
intArgs=(/config_rk_order/), messageType=MPAS_LOG_ERR)
return
endif
endif

! Calculate masks prior to RK loop, but do not update masks within the loop
! to preserve the accuracy of time integration.
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
Comment thread
matthewhoffman marked this conversation as resolved.
err = ior(err, err_tmp)

! *** Start RK loop ***
do rkStage = 1, nRKstages
call mpas_log_write('beginning rk stage $i of $i', &
Expand Down Expand Up @@ -401,9 +405,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + &
(1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:)
thickness = sum(layerThickness, 1)
! Calculate masks after updating thickness
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
err = ior(err, err_tmp)
! Do not calculate masks after updating thickness! We need to keep masks
! constant for now to preserve accuracy of time integration

if (trim(config_thermal_solver) .ne. 'none') then
do iCell = 1, nCells
Expand Down