Skip to content
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

Couple damage to viscosity without modifying stiffness #91

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

Couple damage to viscosity without modifying stiffness #91

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e5ce436
Add damage coupling in velocity solver
trhille Aug 22, 2023
b8a07aa
Fix inadvertent coupling when config_damage_rheology_coupling is false
trhille Aug 24, 2023
ba9c5a9
Don't include damage in uncoupled effectiveViscosity calculation
trhille Aug 25, 2023
9a94b05
Some cleanup
trhille Oct 13, 2023
27e0725
Do not zero out damage for grounded ice
trhille Sep 29, 2023
998e2bd
Disable damage threshold calving for grounded ice
trhille Oct 15, 2023
305eb6a
Write effective stiffness to ascii mesh when coupling damage to visco…
trhille Nov 14, 2023
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
27 changes: 6 additions & 21 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 @@ -3517,11 +3517,10 @@ subroutine li_calculate_damage(domain, err)
! damage is always larger than the Nye value for initialization of damage evolution.

do iCell = 1, nCells
if ((li_mask_is_grounded_ice(cellMask(iCell))) .or. (thickness(iCell) .eq. 0.0_RKIND)) then
if ( thickness(iCell) .eq. 0.0_RKIND ) then
damage(iCell) = 0.0_RKIND
end if
end do
! the damage of grounded ice is kept as 0, as the strain rate calculation is only valid for ice shelf

! Options for initializing damage where ice goes afloat:
if (trim(config_damage_gl_setting) == 'extrapolate') then
Expand Down Expand Up @@ -3554,7 +3553,7 @@ subroutine li_calculate_damage(domain, err)
if (li_mask_is_grounding_line(cellMask(iCell))) then
do iNeighbor = 1, nEdgesOnCell(iCell)
jCell = cellsOnCell(iNeighbor, iCell)
if (li_mask_is_floating_ice(cellMask(jCell))) then
if ( li_mask_is_floating_ice(cellMask(jCell)) .and. (damage(jCell) < damageNye(jCell)) ) then
damage(jCell) = damageNye(jCell)
end if
end do
Expand Down Expand Up @@ -3645,7 +3644,6 @@ subroutine li_finalize_damage_after_advection(domain, err)
type (mpas_pool_type), pointer :: meshPool
type (mpas_pool_type), pointer :: velocityPool
type (mpas_pool_type), pointer :: scratchPool
real(kind=RKIND), pointer :: config_damage_stiffness_min
logical, pointer :: config_damage_rheology_coupling
logical, pointer :: config_preserve_damage
logical, pointer :: config_print_calving_info
Expand All @@ -3667,7 +3665,6 @@ subroutine li_finalize_damage_after_advection(domain, err)
err = 0


call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min)
call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling)
call mpas_pool_get_config(liConfigs, 'config_damage_gl_setting', config_damage_gl_setting)
call mpas_pool_get_config(liConfigs, 'config_preserve_damage', config_preserve_damage)
Expand Down Expand Up @@ -3715,7 +3712,7 @@ subroutine li_finalize_damage_after_advection(domain, err)


do iCell = 1, nCells
if (li_mask_is_grounded_ice(cellMask(iCell)) .or. .not. li_mask_is_ice(cellMask(iCell))) then
if (.not. li_mask_is_ice(cellMask(iCell))) then
damage(iCell) = 0.0_RKIND
end if
end do
Expand Down Expand Up @@ -3751,7 +3748,7 @@ subroutine li_finalize_damage_after_advection(domain, err)
if (li_mask_is_grounding_line(cellMask(iCell))) then
do iNeighbor = 1, nEdgesOnCell(iCell)
jCell = cellsOnCell(iNeighbor, iCell)
if (li_mask_is_floating_ice(cellMask(jCell))) then
if ( li_mask_is_floating_ice(cellMask(jCell)) .and. (damage(jCell)<damageNye(jCell)) ) then
damage(jCell) = damageNye(jCell)
end if
end do
Expand All @@ -3769,18 +3766,6 @@ subroutine li_finalize_damage_after_advection(domain, err)
damage = 1.0_RKIND
end where

if (config_damage_rheology_coupling) then
do iCell = 1, nCells
if (li_mask_is_floating_ice(cellMask(iCell))) then
stiffnessFactor(iCell) = 1.0_RKIND - damage(iCell)
if (stiffnessFactor(iCell) < config_damage_stiffness_min) then
stiffnessFactor(iCell) = config_damage_stiffness_min
end if
end if
end do
end if


block => block % next

enddo
Expand Down Expand Up @@ -3858,11 +3843,11 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d
growMask(:) = 0

! define seed and grow masks for flood fill.
where ( li_mask_is_dynamic_margin(cellMask) .and. (damage .ge. config_damage_calving_threshold) )
where ( li_mask_is_dynamic_margin(cellMask) .and. li_mask_is_floating_ice(cellMask) .and. (damage .ge. config_damage_calving_threshold) )
seedMask = 1
end where

where ( seedMask == 0 .and. (damage .ge. config_damage_calving_threshold) )
where ( seedMask == 0 .and. li_mask_is_floating_ice(cellMask) .and. (damage .ge. config_damage_calving_threshold) )
growMask = 1
end where

Expand Down
21 changes: 18 additions & 3 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -909,12 +909,15 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo
real (kind=RKIND), dimension(:,:), pointer :: flowParamA
real (kind=RKIND), dimension(:), pointer :: effectiveViscosity
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
real (kind=RKIND), dimension(:), pointer :: damage

type (field1dReal), pointer :: meanFlowParamAVar
real (kind=RKIND), dimension(:), pointer :: meanFlowParamA

integer, pointer :: nVertLevels
integer, pointer :: nCells
logical, pointer :: config_damage_rheology_coupling
real (kind=RKIND), pointer :: config_damage_stiffness_min

integer :: iCell

Expand All @@ -926,6 +929,8 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo
call mpas_pool_get_array(velocityPool, 'yvelmean', yvelmean)

call mpas_pool_get_config(liConfigs, 'config_flowLawExponent', config_flowLawExponent)
call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min)
call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling)

call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
Expand All @@ -949,6 +954,7 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo
call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)

call mpas_pool_get_array(geometryPool, 'damage', damage)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'effectiveViscosity', effectiveViscosity)

Expand Down Expand Up @@ -984,9 +990,18 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo
if ( (eEff == 0.0_RKIND) .or. (meanFlowParamA(iCell) == 0.0_RKIND) ) then
effectiveViscosity(iCell) = 0.0_RKIND
else
effectiveViscosity(iCell) = &
0.5_RKIND*stiffnessFactor(iCell)*meanFlowParamA(iCell)**(-1.0_RKIND/config_flowLawExponent)* &
eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent)
if ( config_damage_rheology_coupling ) then
effectiveViscosity(iCell) = &
0.5_RKIND * max(1.0_RKIND - damage(iCell), config_damage_stiffness_min) * &
stiffnessFactor(iCell)*meanFlowParamA(iCell)**(-1.0_RKIND/config_flowLawExponent)* &
eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent)
else
effectiveViscosity(iCell) = &
0.5_RKIND * stiffnessFactor(iCell)*meanFlowParamA(iCell)** &
(-1.0_RKIND/config_flowLawExponent)* &
eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent)
endif

endif
enddo

Expand Down
71 changes: 51 additions & 20 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -449,6 +449,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
normalVelocity, uReconstructX, uReconstructY, uReconstructZ
real (kind=RKIND), dimension(:,:), pointer :: temperature
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
real (kind=RKIND), dimension(:), pointer :: damage
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited
real (kind=RKIND), dimension(:), pointer :: muFriction
Expand All @@ -459,8 +460,10 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
logical, pointer :: config_always_compute_fem_grid
logical, pointer :: config_output_external_velocity_solver_data
logical, pointer :: config_nonconvergence_error
logical, pointer :: config_damage_rheology_coupling
real (kind=RKIND), pointer :: config_ice_density
real (kind=RKIND), pointer :: config_effective_pressure_max
real (kind=RKIND), pointer :: config_damage_stiffness_min
integer, pointer :: anyDynamicVertexMaskChanged
integer, pointer :: dirichletMaskChanged
integer, pointer :: nEdges
Expand All @@ -485,6 +488,8 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
call mpas_pool_get_config(liConfigs, 'config_nonconvergence_error', config_nonconvergence_error)
call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density)
call mpas_pool_get_config(liConfigs, 'config_effective_pressure_max', config_effective_pressure_max)
call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min)
call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling)

! Mesh variables
call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions)
Expand All @@ -502,6 +507,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
call mpas_pool_get_array(geometryPool, 'damage', damage)

! Thermal variables
call mpas_pool_get_array(thermalPool, 'temperature', temperature)
Expand Down Expand Up @@ -595,12 +601,21 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro


call mpas_timer_start("velocity_solver_solve_FO")
call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, &
betaSolve, sfcMassBal, temperature, stiffnessFactor, &
effectivePressureLimited, muFriction, &
uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1
normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values
deltat, albanyVelocityError) ! return values
if ( config_damage_rheology_coupling ) then
call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, &
betaSolve, sfcMassBal, temperature, max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor, &
effectivePressureLimited, muFriction, &
uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1
normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values
deltat, albanyVelocityError) ! return values
else
call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, &
betaSolve, sfcMassBal, temperature, stiffnessFactor, &
effectivePressureLimited, muFriction, &
uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1
normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values
deltat, albanyVelocityError) ! return values
endif
call mpas_timer_stop("velocity_solver_solve_FO")

if (albanyVelocityError == 1) then
Expand Down Expand Up @@ -770,7 +785,8 @@ subroutine li_velocity_external_write_albany_mesh(domain)
!-----------------------------------------------------------------
! local variables
!-----------------------------------------------------------------
logical, pointer :: config_write_albany_ascii_mesh
logical, pointer :: config_write_albany_ascii_mesh, &
config_damage_rheology_coupling
character (len=StrKIND), pointer :: config_velocity_solver
real (kind=RKIND), dimension(:), pointer :: &
bedTopography, lowerSurface, upperSurface, layerThicknessFractions, beta
Expand All @@ -784,12 +800,12 @@ subroutine li_velocity_external_write_albany_mesh(domain)
real (kind=RKIND), dimension(:), pointer :: surfaceAirTemperature, basalHeatFlux
integer, dimension(:), pointer :: vertexMask, cellMask, edgeMask, indexToCellID
real (kind=RKIND), dimension(:,:), pointer :: layerThickness
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor, damage
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: muFriction
integer, dimension(:,:), pointer :: dirichletVelocityMask
type (mpas_pool_type), pointer :: meshPool, geometryPool, thermalPool, observationsPool, velocityPool, scratchPool, hydroPool
real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density
real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density, config_damage_stiffness_min
integer :: iCell
integer :: err

Expand All @@ -810,6 +826,8 @@ subroutine li_velocity_external_write_albany_mesh(domain)
call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density)
call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density)
call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min)
call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling)

if (trim(config_velocity_solver) /= 'FO') then
call mpas_log_write("config_velocity solver needs to be set to 'FO' for config_write_albany_ascii_mesh to work.", &
Expand All @@ -836,6 +854,7 @@ subroutine li_velocity_external_write_albany_mesh(domain)

! Geometry variables
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'damage', damage)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
Expand Down Expand Up @@ -896,17 +915,29 @@ subroutine li_velocity_external_write_albany_mesh(domain)

! call the C++ routine to write the mesh
call mpas_log_write("Writing Albany ASCII mesh.", flushNow=.true.)
call write_ascii_mesh(indexToCellID, bedTopography, lowerSurface, &
beta, temperature, &
surfaceAirTemperature, basalHeatFlux, &
stiffnessFactor, &
effectivePressure, muFriction, &
thickness, thicknessUncertainty, &
sfcMassBal, sfcMassBalUncertainty, &
floatingBasalMassBal, floatingBasalMassBalUncertainty, &
observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, &
observedThicknessTendency, observedThicknessTendencyUncertainty)

if ( config_damage_rheology_coupling ) then
call write_ascii_mesh(indexToCellID, bedTopography, lowerSurface, &
beta, temperature, &
surfaceAirTemperature, basalHeatFlux, &
max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor, &
effectivePressure, muFriction, &
thickness, thicknessUncertainty, &
sfcMassBal, sfcMassBalUncertainty, &
floatingBasalMassBal, floatingBasalMassBalUncertainty, &
observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, &
observedThicknessTendency, observedThicknessTendencyUncertainty)
else
call write_ascii_mesh(indexToCellID, bedTopography, lowerSurface, &
beta, temperature, &
surfaceAirTemperature, basalHeatFlux, &
stiffnessFactor, &
effectivePressure, muFriction, &
thickness, thicknessUncertainty, &
sfcMassBal, sfcMassBalUncertainty, &
floatingBasalMassBal, floatingBasalMassBalUncertainty, &
observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, &
observedThicknessTendency, observedThicknessTendencyUncertainty)
endif
!----

call interface_reset_stdout()
Expand Down