Skip to content

Add embedded Runge-Kutta methods for automatic time-step adjustment #170

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

Add embedded Runge-Kutta methods for automatic time-step adjustment #170

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
9dc0704
Reset masks in RK loop before velocity solve
trhille Aug 29, 2025
dd75c61
Add missing deallocate calls
trhille Mar 4, 2026
b8739a7
Remove outdated comment
trhille Mar 4, 2026
2bdeb34
Add another missing deallocation call
trhille Mar 4, 2026
aff72f2
Control mask updates before velocity solve
trhille Mar 5, 2026
1da7068
Add RK pairs for embedded SSP RK time integration
trhille Mar 12, 2026
11a1141
Only allocate arrays when embedded pairs are enabled
trhille Mar 12, 2026
0faa297
Synchronize embedded norm across ranks
trhille Mar 12, 2026
d1fc644
Do not use damage to calculate error.
trhille Mar 12, 2026
6b83c9e
Make clock advancement consistent with new deltat
trhille Mar 12, 2026
ecae4a4
Fix timekeeping issue in set_timestep
trhille Mar 12, 2026
1177905
Fix deltat update
trhille Mar 12, 2026
bab0fdc
Floor dt before advancing clock
trhille Mar 12, 2026
b822df8
Use global rms error instead of max
trhille Mar 12, 2026
0f341ec
Apply face-melt speed after advection
trhille Mar 12, 2026
064c1b7
Use SSP for embedded pair to RK(4,3)
trhille Mar 13, 2026
f581f11
Change embedded pairs for RK2 and RK(3,3)
trhille Mar 13, 2026
dad6e5b
Remove declaration of unused variable
trhille Mar 13, 2026
9be8aa0
Weight error calculation by cellArea
trhille Mar 13, 2026
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
32 changes: 32 additions & 0 deletions components/mpas-albany-landice/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -558,6 +558,38 @@
description="Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2."
possible_values="3 or 4"
/>
<nml_option name="config_use_embedded_ssprk" type="logical" default_value=".false." units="unitless"
description="Enable embedded SSPRK error estimation and retry-based timestep control for RK3 methods."
possible_values=".true. or .false."
/>
<nml_option name="config_embedded_rk_rtol" type="real" default_value="1.0e-3" units="unitless"
description="Relative tolerance used by embedded SSPRK error control."
possible_values="Any positive real value."
/>
<nml_option name="config_embedded_rk_atol_thickness" type="real" default_value="1.0e-2" units="m"
description="Absolute tolerance for thickness in embedded SSPRK error control."
possible_values="Any positive real value."
/>
<nml_option name="config_embedded_rk_atol_tracer" type="real" default_value="1.0e-4" units="unitless"
description="Absolute tolerance for tracer-like fields in embedded SSPRK error control."
Copy link
Copy Markdown
Author

@trhille trhille Mar 13, 2026

Choose a reason for hiding this comment

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

Yes, it bothered me that co-pilot created an absolute tolerance field that is applied to multiple tracers with different units and magnitudes.

possible_values="Any positive real value."
/>
<nml_option name="config_embedded_rk_safety_factor" type="real" default_value="0.9" units="unitless"
description="Safety factor applied to timestep updates after embedded SSPRK error estimation."
possible_values="Any positive real value less than or equal to 1.0."
/>
<nml_option name="config_embedded_rk_max_growth" type="real" default_value="1.5" units="unitless"
description="Maximum multiplicative timestep increase allowed between embedded SSPRK attempts."
possible_values="Any real value greater than or equal to 1.0."
/>
<nml_option name="config_embedded_rk_min_shrink" type="real" default_value="0.2" units="unitless"
description="Minimum multiplicative timestep shrink factor for rejected embedded SSPRK attempts."
possible_values="Any real value between 0 and 1.0."
/>
<nml_option name="config_embedded_rk_max_retries" type="integer" default_value="3" units="count"
description="Maximum number of embedded SSPRK retry attempts after a rejected step."
possible_values="Any non-negative integer."
/>
<nml_option name="config_adaptive_timestep" type="logical" default_value=".false." units="unitless"
description="Determines if the time step should be adjusted based on the CFL condition or should be steady in time. If true, the config_dt_* options are ignored."
possible_values=".true. or .false."
Expand Down
2 changes: 1 addition & 1 deletion 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 @@ -819,7 +819,7 @@ function li_core_initial_solve(domain) result(err)
call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool)
call mpas_pool_get_array(velocityPool, 'albanyVelocityError', albanyVelocityError)
albanyVelocityError = 0
call li_velocity_solve(domain, solveVelo=solveVelo, err=err_tmp)
call li_velocity_solve(domain, solveVelo=solveVelo, updateMask=.true., err=err_tmp)
err = ior(err, err_tmp)

! Initial calving solution
Expand Down
76 changes: 51 additions & 25 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ module li_iceshelf_melt

!-----------------------------------------------------------------------

subroutine li_face_melt_grounded_ice(domain, err)
subroutine li_face_melt_grounded_ice(domain, err, applyFaceMeltNow, computeFaceMeltNow)
!-----------------------------------------------------------------
! input variables
!-----------------------------------------------------------------
Expand All @@ -88,6 +88,8 @@ subroutine li_face_melt_grounded_ice(domain, err)
!-----------------------------------------------------------------
type (domain_type), intent(inout) :: &
domain !< Input/Output: domain object
logical, intent(in), optional :: applyFaceMeltNow !< If present and false, compute melt fields/CFL only
logical, intent(in), optional :: computeFaceMeltNow !< If present and false, apply existing faceMeltingThickness field

!-----------------------------------------------------------------
! output variables
Expand All @@ -111,9 +113,15 @@ subroutine li_face_melt_grounded_ice(domain, err)
integer, dimension(:), pointer :: cellMask
integer :: err_tmp
logical :: applyToFloating, applyToGrounded, applyToGroundingLine
logical :: doApplyFaceMelt
logical :: doComputeFaceMelt

err = 0
err_tmp = 0
doApplyFaceMelt = .true.
doComputeFaceMelt = .true.
if (present(applyFaceMeltNow)) doApplyFaceMelt = applyFaceMeltNow
if (present(computeFaceMeltNow)) doComputeFaceMelt = computeFaceMeltNow

call mpas_pool_get_config(liConfigs, 'config_front_mass_bal_grounded', config_front_mass_bal_grounded)
call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_float', config_basal_mass_bal_float)
Expand All @@ -133,8 +141,11 @@ subroutine li_face_melt_grounded_ice(domain, err)
if ( trim(config_front_mass_bal_grounded) == 'ismip6' &
.or. trim(config_front_mass_bal_grounded) == 'uniform' ) then
call grounded_face_melt_ismip6(domain, applyToGrounded, &
applyToFloating, applyToGroundingLine, err_tmp)
applyToFloating, applyToGroundingLine, err_tmp, &
computeFaceMeltSpeedNow=doComputeFaceMelt, applyFrontAblationNow=doApplyFaceMelt)
err = ior(err, err_tmp)

if (.not. doApplyFaceMelt) return

Comment on lines 141 to 149
Copy link
Copy Markdown
Author

@trhille trhille Mar 13, 2026

Choose a reason for hiding this comment

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

No, we do not want to convert faceMeltSpeed to faceMeltingThickness until we know the length of the time step, which does not happen until after the RK loop in mpas_li_time_integration_fe_rk.F.

block => domain % blocklist
do while (associated(block))
Expand Down Expand Up @@ -1405,7 +1416,7 @@ end subroutine iceshelf_melt_ismip6


subroutine grounded_face_melt_ismip6(domain, applyToGrounded, &
applyToFloating, applyToGroundingLine, err)
applyToFloating, applyToGroundingLine, err, computeFaceMeltSpeedNow, applyFrontAblationNow)

use li_calving
use li_constants, only: oceanFreezingTempDepthDependence
Expand All @@ -1416,6 +1427,8 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, &
type (domain_type), intent(inout) :: &
domain !< Input/Output: domain object
logical, intent(in) :: applyToFloating, applyToGrounded, applyToGroundingLine
logical, intent(in), optional :: computeFaceMeltSpeedNow
logical, intent(in), optional :: applyFrontAblationNow
!-----------------------------------------------------------------
! output variables
!-----------------------------------------------------------------
Expand Down Expand Up @@ -1464,11 +1477,16 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, &
logical, pointer :: config_calculate_thermal_forcing
logical, pointer :: config_ocean_data_extrapolation
integer, pointer :: config_ocean_data_extrap_ncells_extra
logical :: doComputeFaceMeltSpeed, doApplyFrontAblation
real (kind=RKIND) :: dzAccumulated, dz
integer :: err_tmp, kk, iLayer

err = 0
call mpas_log_write('Starting face melt routine')
doComputeFaceMeltSpeed = .true.
doApplyFrontAblation = .true.
if (present(computeFaceMeltSpeedNow)) doComputeFaceMeltSpeed = computeFaceMeltSpeedNow
if (present(applyFrontAblationNow)) doApplyFrontAblation = applyFrontAblationNow

! Get sea level, bedTopography, ice density
call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
Expand Down Expand Up @@ -1511,7 +1529,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, &
call mpas_pool_get_array(geometryPool, 'faceMeltingCFLdt', faceMeltingCFLdt)
call mpas_pool_get_array(geometryPool, 'dtFaceMeltingCFLratio', dtFaceMeltingCFLratio)

if ( config_use_3d_thermal_forcing_for_face_melt ) then
if (doComputeFaceMeltSpeed .and. config_use_3d_thermal_forcing_for_face_melt ) then
call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers)
call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing)
call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', zOcean)
Expand Down Expand Up @@ -1624,29 +1642,37 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, &
TFocean = TFocean + ismip6shelfMelt_deltaT !apply dT correction
endif

faceMeltSpeed(:) = 0.0_RKIND
allocate(faceMeltSpeedVertAvg(nCells+1))
faceMeltSpeedVertAvg(:) = 0.0_RKIND
if (doComputeFaceMeltSpeed) then
faceMeltSpeed(:) = 0.0_RKIND

! Calculate face melting for each cell
do iCell = 1, nCells
! Calculate face-melt speed for each cell
do iCell = 1, nCells

if ( bedTopography(iCell) < 0.0_RKIND ) then
waterDepth = seaLevel - bedTopography(iCell)
else
waterDepth = 0.0_RKIND
endif
if ( bedTopography(iCell) < 0.0_RKIND ) then
waterDepth = seaLevel - bedTopography(iCell)
else
waterDepth = 0.0_RKIND
endif

if (trim(config_front_mass_bal_grounded) == 'ismip6') then
! Calculate ice front melt rate at each cell
faceMeltSpeed(iCell) = (aSubglacial * waterDepth * & ! m s^-1
(ismip6Runoff(iCell) / rhow * secPerDay)**alphaSubglacial &
+ B) * max(0.0_RKIND, TFocean(iCell) + addTFocean)**betaTF / secPerDay
elseif (trim(config_front_mass_bal_grounded) == 'uniform') then
faceMeltSpeed(iCell) = config_uniform_face_melt_rate
endif
if (trim(config_front_mass_bal_grounded) == 'ismip6') then
! Calculate ice front melt rate at each cell
faceMeltSpeed(iCell) = (aSubglacial * waterDepth * & ! m s^-1
(ismip6Runoff(iCell) / rhow * secPerDay)**alphaSubglacial &
+ B) * max(0.0_RKIND, TFocean(iCell) + addTFocean)**betaTF / secPerDay
elseif (trim(config_front_mass_bal_grounded) == 'uniform') then
faceMeltSpeed(iCell) = config_uniform_face_melt_rate
endif

enddo
enddo
endif

if (.not. doApplyFrontAblation) then
block => block % next
cycle
endif

Comment on lines +1669 to +1673
Copy link
Copy Markdown
Author

@trhille trhille Mar 13, 2026

Choose a reason for hiding this comment

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

The CFL calculation might be a legitimate concern, but we definitely do not want to calculate faceMeltingThickness yet.

allocate(faceMeltSpeedVertAvg(nCells+1))
faceMeltSpeedVertAvg(:) = 0.0_RKIND

where ( (thickness > 0.0_RKIND) .and. (lowerSurface < 0.0_RKIND) )
faceMeltSpeedVertAvg(:) = faceMeltSpeed(:) * abs(lowerSurface(:)) / thickness(:)
Expand All @@ -1660,12 +1686,12 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, &
maxDt=faceMeltingCFLdt, CFLratio=dtFaceMeltingCFLratio, err=err_tmp)
err = ior(err, err_tmp)

deallocate(faceMeltSpeedVertAvg)

Comment on lines 1674 to +1690

block => block % next
enddo ! associated(block)

deallocate(faceMeltSpeedVertAvg)

end subroutine grounded_face_melt_ismip6
!-----------------------------------------------------------------------

Expand Down
Loading