Skip to content

Commit

Permalink
Check master dt/SGH coupling interval compatibility
Browse files Browse the repository at this point in the history 
Introduces a new subroutine to the SGH model that ensures
config_adaptive_interval_force_interval divides evenly into
config_SGH_coupling_interval
  • Loading branch information
@alexolinhager
alexolinhager committed Sep 10, 2024
1 parent 5ea7382 commit 785649e
Show file tree
Hide file tree
Showing 2 changed files with 105 additions and 8 deletions.
4 changes: 3 additions & 1 deletion components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,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
109 changes: 102 additions & 7 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,9 +117,12 @@ subroutine li_SGH_init(domain, err)
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 @@ -133,10 +136,17 @@ subroutine li_SGH_init(domain, err)
return
endif

! check SGH coupling interval
! Check SGH coupling interval is compatible with master timestep and restart 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, timeStep=config_SGH_coupling_interval, ierr=err_tmp)
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)
Expand All @@ -145,7 +155,7 @@ subroutine li_SGH_init(domain, err)
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(domainb%clock, 'sghCouplingInterval', ierr=err_tmp)
call mpas_reset_clock_alarm(domain%clock, 'sghCouplingInterval', ierr=err_tmp)
err = ior(err, err_tmp)
endif

Expand Down Expand Up @@ -333,7 +343,7 @@ subroutine li_SGH_solve(domain, err)
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 :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureTimeCycleAvg
real (kind=RKIND), dimension(:), allocatable :: cellJunk
integer, dimension(:), pointer :: nEdgesOnCell
Expand Down Expand Up @@ -747,7 +757,7 @@ subroutine li_SGH_solve(domain, err)
! set up coupling interval
if (mpas_is_alarm_ringing(domain % clock, 'sghCouplingInterval', ierr=err_tmp)) then
err = ior(err, err_tmp)
block => doimain % blocklist
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
Expand All @@ -761,7 +771,7 @@ subroutine li_SGH_solve(domain, err)
err = ior(err, err_tmp)
endif
call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/)
call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/))
! === error check
if (err > 0) then
Expand Down Expand Up @@ -2776,4 +2786,89 @@ 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.
!
!-----------------------------------------------------------------------
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_set_timeInterval(coupling_interval, timeString=config_SGH_coupling_interval, ierr=err_tmp)
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

0 comments on commit 785649e

Please sign in to comment.