Skip to content

Set SGH coupling interval #125

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
alexolinhager wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Set SGH coupling interval #125

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1247e82
Establish initial alarm framework
alexolinhager Sep 3, 2024
5ea7382
Complete effectivePressure averaging over coupling interval
alexolinhager Sep 4, 2024
727dcd9
Check master dt/SGH coupling interval compatibility
alexolinhager Sep 5, 2024
f123500
Remove effectivePressure set to zero at timeinterval
alexolinhager Sep 16, 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
16 changes: 12 additions & 4 deletions components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,11 @@
description="The maximum allowable time step in seconds. If the allowable time step determined by the adaptive CFL calculation is longer than this, then the model will specify config_SGH_max_adaptive_timestep as the time step instead. Defaults to 100 years (in seconds)."
possible_values="Any non-negative real value."
/>
<nml_option name="config_SGH_tangent_slope_calculation" type="character" default_value="from_normal_slope" units="unitless"
<nml_option name="config_SGH_coupling_interval" type="character" default_value="0000-00-00_00:00:01" units="unitless"
description="Time interval at which the SGH model is called by MALI. If set to an interval less than the master timestep, config_SGH_coupling_interval will be superceded by the master timestep. If greater than the master timestep, config_adaptive_timestep_force_interval will be set equal to config_SGH_coupling_interval."
Copy link
Copy Markdown

@matthewhoffman matthewhoffman Sep 14, 2024

Choose a reason for hiding this comment

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

Rather than modifying config_adaptive_timestep_force_interval, it might be preferable to throw an error if it's incompatible with the new option. I'm thinking this because config_adaptive_timestep_force_interval is used for other things too and it might get wonky if some parts of the code or adjusting it.

possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS'"
/>
<nml_option name="config_SGH_tangent_slope_calculation" type="character" default_value="from_normal_slope" units="unitless"
description="Selection of the method for calculating the tangent component of slope at edges.
'from_vertex_barycentric' interpolates scalar values from cell centers to vertices using the barycentric interpolation routine in operators (mpas_cells_to_points_using_baryweights) and then calculates the slope between vertices. It works for obtuse triangles, but will not work correctly across the edges of periodic meshes.
'from_vertex_barycentric_kiteareas' interpolates scalar values from cell centers to vertices using barycentric interpolation based on kiterea values and then calculates the slope between vertices. It will work across the edges of periodic meshes, but will not work correctly for obtuse triangles.
Expand Down Expand Up @@ -180,8 +184,10 @@
<!-- convenience variables -->
<var name="effectivePressureSGH" type="real" dimensions="nCells Time" units="Pa"
description="effective ice pressure in subglacial hydrology system" />
<var name="effectivePressure" type="real" dimensions="nCells Time" units="Pa" packages="higherOrderVelocity"
description="Effective pressure used in basal friction calculation. If subglacial hydrology model is active, this will be effectivePressureSGH averaged over the subglacial hydrology model timestepping subcycles. If subglacial hydrology model is inactive, this will come from a file or a parameterization."/>
<var name="effectivePressure" type="real" dimensions="nCells Time" units="Pa" packages="higherOrderVelocity" default_value="-1e36"
description="Effective pressure used in basal friction calculation. If subglacial hydrology model is active, this will be effectivePressureSGH averaged over the interval defined by config_SGH_coupling_interval. If subglacial hydrology model is inactive, this will come from a file or a parameterization."/>
<var name="effectivePressureTimeCycleAvg" type="real" dimensions="nCells Time" units="Pa"
description="effectivePressureSGH averaged over the subglacial hydrology model timestepping subcycle." />
<var name="hydropotential" type="real" dimensions="nCells Time" units="Pa"
description="hydropotential in subglacial hydrology system" />
<var name="waterFlux" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
Expand Down Expand Up @@ -246,7 +252,9 @@
description="time step length limited by pressure equation scheme in subglacial hydrology system" />
<var name="deltatSGH" type="real" dimensions="Time" units="s"
description="time step used for evolving subglacial hydrology system" />
<var name="distGroundingLineDischargeCell" type="real" dimensions="Time nCells" units="m^{3} s^{-1}"
<var name="deltatSGHaccumulated" type="real" dimensions="Time" units="s"
description="accumulated deltatSGH throughout SGH coupling interval" />
<var name="distGroundingLineDischargeCell" type="real" dimensions="Time nCells" units="m^{3} s^{-1}"
description="distributed discharge across the grounding line, summed from grounding line edges to adjacent ungrounded cell. Values from all edges are summed if multiple grounding line edges border a single ungrounded cell" />
<var name="totalGroundingLineDischargeCell" type="real" dimensions="Time nCells" units="m^{3} s^{-1}"
description="total (channel + dist.) discharge across the grounding line, summed from grounding line edges to adjacent ungrounded cell. Values from all edges are summed if multiple grounding line edges border a single ungrounded cell" />
Expand Down
181 changes: 165 additions & 16 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,15 +111,20 @@ subroutine li_SGH_init(domain, err)
real (kind=RKIND), dimension(:), pointer :: bedTopography
real (kind=RKIND), dimension(:), pointer :: iceThicknessHydro
real (kind=RKIND), dimension(:), pointer :: hydropotential
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND),dimension(:), pointer :: effectivePressure
integer, dimension(:), pointer :: cellMask
real (kind=RKIND), pointer :: tillMax
real (kind=RKIND), pointer :: rhoi, rhoo
logical, pointer :: config_do_restart
real (kind=RKIND), pointer :: config_sea_level
integer, pointer :: config_num_halos
character (len=StrKIND), pointer :: config_SGH_coupling_interval
type (MPAS_TimeInterval_type) :: sgh_coupling_interval
character (len=StrKIND), pointer :: simulationStartTime
type (MPAS_Time_type) :: simulationStartTime_timeType
integer :: err_tmp


err = 0
err_tmp = 0

Expand All @@ -132,8 +137,31 @@ subroutine li_SGH_init(domain, err)
! If SGH is not active, skip everything
return
endif



! Check SGH coupling interval is compatible with master timestep and restart interval
Copy link
Copy Markdown

@matthewhoffman matthewhoffman Sep 14, 2024

Choose a reason for hiding this comment

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

Suggested change
! Check SGH coupling interval is compatible with master timestep and restart interval
! Check SGH coupling interval is compatible with master timestep and timestepping force interval

call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
call check_SGH_coupling_interval(meshPool, err_tmp)
err = ior(err, err_tmp)
if (err /= 0) then
call mpas_log_write("Error occurred in check_SGH_coupling_interval. ", MPAS_LOG_ERR)
endif

! Initiate coupling interval timer
call mpas_pool_get_config(liConfigs, 'config_SGH_coupling_interval', config_SGH_coupling_interval)
call mpas_set_timeInterval(sgh_coupling_interval, timeString=config_SGH_coupling_interval, ierr=err_tmp)
err = ior(err, err_tmp)
call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime)
call mpas_set_time(simulationStartTime_timeType, dateTimeString=simulationStartTime, ierr=err_tmp)
err = ior(err, err_tmp)
call mpas_add_clock_alarm(domain%clock,'sghCouplingInterval', alarmTime=simulationStartTime_timeType, alarmTimeInterval=sgh_coupling_interval, ierr=err_tmp)
err = ior(err, err_tmp)
if (mpas_is_alarm_ringing(domain%clock, 'sghCouplingInterval', ierr=err_tmp)) then
err = ior(err, err_tmp)
call mpas_reset_clock_alarm(domain%clock, 'sghCouplingInterval', ierr=err_tmp)
err = ior(err, err_tmp)
endif

! begin hydro init
call mpas_timer_start("hydro init")

call mpas_log_write('Beginning subglacial hydro init.')
Expand Down Expand Up @@ -239,7 +267,22 @@ subroutine li_SGH_init(domain, err)
call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'waterPressureSmooth')
call mpas_timer_stop("halo updates")


! initialize effectivePressure
block => domain % blocklist
do while (associated(block))

call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)

if (all(effectivePressure == -1e36_RKIND)) then
effectivePressure = effectivePressureSGH
endif

block => block % next
end do

! === error check
if (err > 0) then
call mpas_log_write("An error has occurred in li_SGH_init.", MPAS_LOG_ERR)
Expand Down Expand Up @@ -316,23 +359,24 @@ subroutine li_SGH_solve(domain, err)
real (kind=RKIND), dimension(:), pointer :: waterVelocity
real (kind=RKIND), dimension(:), pointer :: waterVelocityCellX
real (kind=RKIND), dimension(:), pointer :: waterVelocityCellY
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: effectivePressureTimeCycleAvg
real (kind=RKIND), dimension(:), allocatable :: cellJunk
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge
integer, dimension(:,:), pointer :: edgeSignOnCell
integer, dimension(:), pointer :: cellMask
real (kind=RKIND), pointer :: deltatSGH
real (kind=RKIND), pointer :: deltatSGHaccumulated
real (kind=RKIND), pointer :: masterDeltat
real (kind=RKIND), pointer :: Cd
real (kind=RKIND), pointer :: tillMax
integer, pointer :: nCellsSolve
integer, pointer :: nCells
integer :: iCell, iEdge, iEdgeOnCell
real (kind=RKIND) :: timeLeft ! subcycling time remaining until master dt is reached
real (kind=RKIND) :: deltatSGHaccumulated
integer :: numSubCycles ! number of subcycles
integer :: err_tmp

Expand Down Expand Up @@ -380,25 +424,27 @@ subroutine li_SGH_solve(domain, err)
call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'iceThicknessHydro')
call mpas_timer_stop("halo updates")

! initialize while loop
call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) ! can get from any block
call mpas_pool_get_array(meshPool, 'deltat', masterDeltat)
timeLeft = masterDeltat ! in seconds
numSubCycles = 0

! initialize accumulated time averages/sums throughout timecycle
block => domain % blocklist
do while (associated(block))

call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
effectivePressure = 0.0_RKIND
call mpas_pool_get_array(hydroPool, 'effectivePressureTimeCycleAvg', effectivePressureTimeCycleAvg)
call mpas_pool_get_array(hydroPool, 'deltatSGHaccumulated', deltatSGHaccumulated)

effectivePressureTimeCycleAvg = 0.0_RKIND
deltatSGHaccumulated = 0.0_RKIND

block => block % next
end do


enddo

! =============
! =============
! =============
Expand Down Expand Up @@ -704,10 +750,11 @@ subroutine li_SGH_solve(domain, err)

call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'effectivePressureTimeCycleAvg', effectivePressureTimeCycleAvg)
call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
call mpas_pool_get_array(hydroPool, 'deltatSGHaccumulated', deltatSGHaccumulated)

effectivePressure = (effectivePressure * deltatSGHaccumulated + effectivePressureSGH * deltatSGH) / (deltatSGHaccumulated + deltatSGH)
effectivePressureTimeCycleAvg = (effectivePressureTimeCycleAvg * deltatSGHaccumulated + effectivePressureSGH * deltatSGH) / (deltatSGHaccumulated + deltatSGH)

block => block % next
end do
Expand All @@ -722,8 +769,24 @@ subroutine li_SGH_solve(domain, err)
! =============
! =============

call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/))
! set up coupling interval
if (mpas_is_alarm_ringing(domain % clock, 'sghCouplingInterval', ierr=err_tmp)) then
err = ior(err, err_tmp)
block => domain % blocklist
do while (associated(block))

call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'effectivePressureTimeCycleAvg', effectivePressureTimeCycleAvg)

effectivePressure = effectivePressureTimeCycleAvg
block => block % next
end do
call mpas_reset_clock_alarm(domain % clock, 'sghCouplingInterval', ierr=err_tmp)
err = ior(err, err_tmp)
endif

call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/))

! === error check
if (err > 0) then
Expand Down Expand Up @@ -2738,4 +2801,90 @@ subroutine calc_waterPressureSmooth(block,err)
deallocate(pressureField)
end subroutine calc_waterPressureSmooth

!***********************************************************************
!
! routine check_SGH_coupling_interval
!
!> \brief Perform various checks on the SGH coupling interval setting
!> \author Matt Hoffman, modified by Alex Hager
!> \date Sep 2024
!> \details
!> This routine checks that the SGH coupling interval is an even multiple
!> of the adaptive timestep force inverval and divides evenly into the
!> restart interval.
Copy link
Copy Markdown

@matthewhoffman matthewhoffman Sep 14, 2024

Choose a reason for hiding this comment

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

I think we settled on not checking the restart interval, because we don't need the SGH coupling interval to divide evenly into the restart interval if we add the effective pressure accumulator fields as restart fields. I realize this comment might be a holdover from where this was copied from, so just flagging this for follow up.

!
!-----------------------------------------------------------------------
subroutine check_SGH_coupling_interval(meshPool, err)

use mpas_timekeeping
use mpas_stream_manager
use mpas_derived_types, only : MPAS_STREAM_PROPERTY_RECORD_INTV

type (mpas_pool_type), intent(in) :: meshPool !< mesh information
integer, intent(out) :: err

! local variables
character (len=StrKIND), pointer :: config_SGH_coupling_interval
logical, pointer :: config_adaptive_timestep
character (len=StrKIND), pointer :: config_adaptive_timestep_force_interval, config_dt
type (MPAS_TimeInterval_Type) :: coupling_interval, force_interval, dt_interval, zero_interval
type (MPAS_Time_Type) :: start_time
character (len=StrKIND), pointer :: simulationStartTime
integer (kind=I8KIND) :: n_intervals
type (MPAS_TimeInterval_Type) :: remainder
integer :: err_tmp

err = 0
err_tmp = 0

call mpas_log_write("")
call mpas_log_write("-- Checking consistency of config_SGH_coupling_interval and other settings --")

! define zero interval for comparing against below
call mpas_set_timeInterval(zero_interval, dt = 0.0_RKIND)
! get start time as a reference time
call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime)
call mpas_set_time(start_time, dateTimeString=simulationStartTime)
! define SGH coupling interval as a timeInterval type
call mpas_pool_get_config(liConfigs, 'config_SGH_coupling_interval', config_SGH_coupling_interval)
call mpas_set_timeInterval(coupling_interval, timeString=config_SGH_coupling_interval, ierr=err_tmp)
Copy link
Copy Markdown

@matthewhoffman matthewhoffman Sep 14, 2024

Choose a reason for hiding this comment

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

need to retrieve config_SGH_coupling_interval before using it

err = ior(err, err_tmp)

call mpas_pool_get_config(liConfigs, "config_adaptive_timestep", config_adaptive_timestep)
if (config_adaptive_timestep) then
! for adaptive dt, check that config_adaptive_timestep_force_interval divides evenly into config_SGH_coupling_interval
call mpas_pool_get_config(liConfigs, "config_adaptive_timestep_force_interval", config_adaptive_timestep_force_interval)
call mpas_set_timeInterval(force_interval, timeString=config_adaptive_timestep_force_interval, ierr=err_tmp)
err = ior(err, err_tmp)
call mpas_interval_division(start_time, coupling_interval, force_interval, n_intervals, remainder)
if (remainder .EQ. zero_interval) then
call mpas_log_write("config_adaptive_timestep_force_interval divides into config_SGH_coupling_interval $i times " // &
"with no remainder - check passes", intArgs=(/int(n_intervals)/))
else
call mpas_log_write("config_adaptive_timestep_force_interval divides into config_SGH_coupling_interval $i times " // &
"with nonzero remainder", MPAS_LOG_ERR, intArgs=(/int(n_intervals)/))
err = ior(err, 1)
endif
else
! For fixed dt, check that dt divides evenly into config_SGH_coupling_interval
call mpas_pool_get_config(liConfigs, "config_dt", config_dt)
call mpas_set_timeInterval(dt_interval, timeString=config_dt, ierr=err_tmp)
err = ior(err, err_tmp)
call mpas_interval_division(start_time, coupling_interval, dt_interval, n_intervals, remainder)
if (remainder .EQ. zero_interval) then
call mpas_log_write("config_dt divides into config_SGH_coupling_interval $i times " // &
"with no remainder - check passes", intArgs=(/int(n_intervals)/))
else
call mpas_log_write("config_dt divides into config_SGH_coupling_interval $i times " // &
"with nonzero remainder", MPAS_LOG_ERR, intArgs=(/int(n_intervals)/))
err = ior(err, 1)
endif
endif

call mpas_log_write("")

!--------------------------------------------------------------------
end subroutine check_SGH_coupling_interval


end module li_subglacial_hydro