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

Modify stiffness using damage change #55

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

Modify stiffness using damage change #55

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
8a091ed
Initialize damage using stiffnessFactor
trhille Sep 23, 2022
1e1880e
Fix call to mpas_pool_get_subpool
trhille Sep 23, 2022
f30ed87
Modify stiffnessFactor using change in damage
trhille Sep 23, 2022
7388df0
Update Registry for new damage options
trhille Sep 23, 2022
1a478a0
Do not initialize damage from stiffnessFactor on restarts
trhille Dec 15, 2022
9a567de
Update halos for damage before applying to stiffnessFactor
trhille Dec 16, 2022
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
6 changes: 5 additions & 1 deletion components/mpas-albany-landice/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,7 +268,11 @@
description="If true, then the value of 'stiffnessFactor' is coupled to damage evolution, i.e. 'stiffnessFactor' = 1-damage."
possible_values=".true. or .false."
/>
<nml_option name="config_damage_gl_setting" type="character" default_value="nye" units="unitless"
<nml_option name="config_initialize_damage_with_stiffnessFactor" type="logical" default_value=".false." units="unitless"
description="If true, then initialize damage to max((1 - stiffnessFactor), 0)"
possible_values=".true. or .false."
/>
<nml_option name="config_damage_gl_setting" type="character" default_value="nye" units="unitless"
description="Selection of the method for initializing damage in the first floating grid cells past the grounding line: 'nye' uses Nye's zero stress criteria (Nye, 1957, Proc. R. Soc. A, 239) for defining the damage value, 'extrapolate' extrapolates values from the floating ice downstream back to the first floating cells."
possible_values="'nye', 'extrapolate'"
/>
Expand Down
20 changes: 18 additions & 2 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 @@ -3535,6 +3535,11 @@ subroutine li_calculate_damage(domain, err)
err = ior(err, 1)
endif

call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'damage')
call mpas_dmpar_field_halo_exch(domain, 'damageNye')
call mpas_timer_stop("halo updates")

if (config_print_calving_info) then
! End of routine accounting/reporting
localMinInfo(1) = minval(damageSource)
Expand Down Expand Up @@ -3616,11 +3621,13 @@ subroutine li_finalize_damage_after_advection(domain, err)
type (mpas_pool_type), pointer :: velocityPool
type (mpas_pool_type), pointer :: scratchPool
real(kind=RKIND), pointer :: config_damage_stiffness_min
real(kind=RKIND), pointer :: deltat
logical, pointer :: config_damage_rheology_coupling
logical, pointer :: config_preserve_damage
logical, pointer :: config_print_calving_info
character (len=StrKIND), pointer :: config_damage_gl_setting
real (kind=RKIND), dimension(:), pointer :: damage
real (kind=RKIND), dimension(:), pointer :: ddamagedt
real (kind=RKIND), dimension(:), pointer :: damageMax
real (kind=RKIND), dimension(:), pointer :: damageNye
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
Expand Down Expand Up @@ -3653,11 +3660,12 @@ subroutine li_finalize_damage_after_advection(domain, err)

! get fields
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)

call mpas_pool_get_array(meshPool, 'deltat', deltat)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'damage', damage)
call mpas_pool_get_array(geometryPool, 'ddamagedt', ddamagedt)
call mpas_pool_get_array(geometryPool, 'damageMax', damageMax)
call mpas_pool_get_array(geometryPool, 'damageNye', damageNye)
call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
Expand Down Expand Up @@ -3739,10 +3747,18 @@ subroutine li_finalize_damage_after_advection(domain, err)
damage = 1.0_RKIND
end where

! Halo update for damage before applying to stiffnessFactor
! NOTE: THIS WILL NOT WORK ON MULTIPLE BLOCKS PER PROCESSOR
call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'damage')
call mpas_dmpar_field_halo_exch(domain, 'damageNye')
call mpas_timer_stop("halo updates")

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)
stiffnessFactor(iCell) = stiffnessFactor(iCell) - &
(ddamagedt(iCell) * deltat)
if (stiffnessFactor(iCell) < config_damage_stiffness_min) then
stiffnessFactor(iCell) = config_damage_stiffness_min
end if
Expand Down
17 changes: 17 additions & 0 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_core.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -749,6 +749,7 @@ function li_core_initial_solve(domain) result(err)
type (mpas_pool_type), pointer :: geometryPool, meshPool, velocityPool
logical, pointer :: config_do_restart, config_write_output_on_startup, config_write_stats_on_startup
character(len=StrKIND), pointer :: config_velocity_solver
logical, pointer :: config_initialize_damage_with_stiffnessFactor
! Variables needed for printing timestamps
type (MPAS_Time_Type) :: currTime
character(len=StrKIND) :: timeStamp
Expand All @@ -759,6 +760,8 @@ function li_core_initial_solve(domain) result(err)
logical :: solveVelo

integer, dimension(:), pointer :: vertexMask
real (kind=RKIND), dimension(:), pointer :: damage
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor


err = 0
Expand All @@ -770,14 +773,28 @@ function li_core_initial_solve(domain) result(err)
call mpas_pool_get_config(liConfigs, 'config_write_output_on_startup', config_write_output_on_startup)
call mpas_pool_get_config(liConfigs, 'config_write_stats_on_startup', config_write_stats_on_startup)
call mpas_pool_get_config(liConfigs, 'config_velocity_solver', config_velocity_solver)
call mpas_pool_get_config(liConfigs, 'config_initialize_damage_with_stiffnessFactor', config_initialize_damage_with_stiffnessFactor)
call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool)

currTime = mpas_get_clock_time(domain % clock, MPAS_NOW, err_tmp)
err = ior(err, err_tmp)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=err_tmp)
err = ior(err, err_tmp)
call mpas_log_write('Initial timestep ' // trim(timeStamp))

! initialize damage using stiffnessFactor
call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
call mpas_pool_get_array(geometryPool, 'damage', damage)
if ( config_initialize_damage_with_stiffnessFactor .and. (.not. config_do_restart) ) then
damage(:) = 1.0_RKIND - stiffnessFactor(:)

where ( damage < 0.0_RKIND )
damage = 0.0_RKIND
end where

end if
! ===
! === Calculate Initial state
! ===
Expand Down