From 8a091edc052f04a005ea777f8ee427dc5a9de610 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 23 Sep 2022 09:05:23 -0600 Subject: [PATCH 1/6] Initialize damage using stiffnessFactor Previously with damage-rheology coupling there was no stiffnessFactor constraint on initial damage. Now use damage = (1 - stiffnessFactor) to initialize. --- .../src/mode_forward/mpas_li_core.F | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index 10ee1e94b305..947fae6718aa 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -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 @@ -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 @@ -770,6 +773,10 @@ 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(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) + call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) currTime = mpas_get_clock_time(domain % clock, MPAS_NOW, err_tmp) err = ior(err, err_tmp) @@ -777,7 +784,17 @@ function li_core_initial_solve(domain) result(err) 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 ) then + damage(:) = 1.0_RKIND - stiffnessFactor(:) + where ( damage < 0.0_RKIND ) + damage = 0.0_RKIND + end where + + end if ! === ! === Calculate Initial state ! === From 1e1880e33bf41a81a93d3871c23e84fd3bdefeaf Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 23 Sep 2022 11:48:55 -0600 Subject: [PATCH 2/6] Fix call to mpas_pool_get_subpool Need to use domain % blocklist % structs instead of block % structs when outside loop over blocks. --- .../mpas-albany-landice/src/mode_forward/mpas_li_core.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index 947fae6718aa..55c4db4107e4 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -774,9 +774,9 @@ function li_core_initial_solve(domain) result(err) 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(block % structs, 'mesh', meshPool) - call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) + 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) From f30ed87367d857e2bedb89f56a85c3e9e045e4ce Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 23 Sep 2022 12:48:27 -0600 Subject: [PATCH 3/6] Modify stiffnessFactor using change in damage Because stiffnessFactor is often >1 from the optimization, we need to update it using change in damage, rather than setting it to 1 - damage, which previously always resulted in stiffness <1. --- .../src/mode_forward/mpas_li_calving.F | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index d6fee67ee637..b90dcc3901b6 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -3616,11 +3616,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 @@ -3653,11 +3655,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) @@ -3742,7 +3745,8 @@ subroutine li_finalize_damage_after_advection(domain, err) 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 From 7388df02751fa350d5a30348f0df36c3636c82c1 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 23 Sep 2022 14:24:17 -0600 Subject: [PATCH 4/6] Update Registry for new damage options Forgot to add Registry to previous commits! --- components/mpas-albany-landice/src/Registry.xml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index a6cf3e947a81..3ab91cc321a6 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -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." /> - + From 1a478a06baae20dffd67a86172c61e8916de202c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 15 Dec 2022 14:36:01 -0800 Subject: [PATCH 5/6] Do not initialize damage from stiffnessFactor on restarts On a restart, just use the damage field in the restart file, rather than calculating it from stiffnessFactor. --- components/mpas-albany-landice/src/mode_forward/mpas_li_core.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index 55c4db4107e4..ca9069419764 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -787,7 +787,7 @@ function li_core_initial_solve(domain) result(err) ! 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 ) then + if ( config_initialize_damage_with_stiffnessFactor .and. (.not. config_do_restart) ) then damage(:) = 1.0_RKIND - stiffnessFactor(:) where ( damage < 0.0_RKIND ) From 9a567ded918df63a3a1195585961d97836173d0c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 16 Dec 2022 10:20:45 -0800 Subject: [PATCH 6/6] Update halos for damage before applying to stiffnessFactor Update halos for damage before applying to stiffnessFactor. Testing showed imprinting of the decomposition on damage and stiffnessFactor fields before applying this halo update. --- .../src/mode_forward/mpas_li_calving.F | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index b90dcc3901b6..f63bf787b0c5 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -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) @@ -3742,6 +3747,13 @@ 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