Skip to content

Commit

Permalink
Add option for 4-stage SSPRK3
Browse files Browse the repository at this point in the history 
Add option for 4-stage SSPRK3, which involves four velocities solves
instead of three, but also allows for a CFL fraction of 2. To use
the 4-stage formulation, set 'config_rk_order = 3' and
'config_rk3_stages = 4'.
  • Loading branch information
@trhille
trhille committed Sep 20, 2023
1 parent 1f6574c commit 4f69cfe
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 15 deletions.
4 changes: 4 additions & 0 deletions components/mpas-albany-landice/src/Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -492,6 +492,10 @@
description="Order of Runge-Kutta time integration to use. A value of 1 is equivalent to forward euler. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods."
possible_values="1, 2, 3"
/>
<nml_option name="config_rk3_stages" type="integer" default_value="3" units="stages"
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_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
58 changes: 43 additions & 15 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
character (len=StrKIND), pointer :: config_thickness_advection
character (len=StrKIND), pointer :: config_thermal_solver
character (len=StrKIND), pointer :: config_time_integration
integer, pointer :: config_rk_order
integer, pointer :: config_rk_order, config_rk3_stages
integer :: rkStage, iCell, iTracer, k
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, &
temperature, &
Expand All @@ -133,7 +133,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
real (kind=RKIND), pointer :: deltat
real (kind=RKIND) :: deltatFull
real (kind=RKIND), dimension(4) :: rkSubstepWeights
real (kind=RKIND), dimension(3) :: rkSSPweights
real (kind=RKIND), dimension(4) :: rkSSPweights
integer :: nRKstages

err = 0
Expand All @@ -145,6 +145,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection)
call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver)
call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order)
call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages)
call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration)

call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
Expand Down Expand Up @@ -262,16 +263,33 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
rkSSPweights(1) = 1.0_RKIND
rkSSPweights(2) = 0.5_RKIND
rkSSPweights(3) = 0.0_RKIND
rkSSPweights(4) = 0.0_RKIND
elseif ( (trim(config_time_integration) == 'runge_kutta') .and. &
(config_rk_order == 3) ) then
! use Strong Stability Preserving RK3
nRKstages = 3
rkSubstepWeights(:) = 1.0_RKIND
rkSSPweights(1) = 1.0_RKIND
rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND
rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND
! Add option for 4-stage SSP-RK3? Allows for CFL=2.
if (config_rk3_stages == 3) then
! use Strong Stability Preserving RK3
nRKstages = 3
rkSubstepWeights(:) = 1.0_RKIND
rkSSPweights(1) = 1.0_RKIND
rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND
rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND
rkSSPweights(4) = 0.0_RKIND
! 4-stage SSP-RK3? Allows for CFL=2.
elseif (config_rk3_stages == 4) then
nRKstages = 4
rkSubstepWeights(:) = 0.5_RKIND
rkSSPweights(1) = 1.0_RKIND
rkSSPweights(2) = 0.0_RKIND
rkSSPweights(3) = 2.0_RKIND / 3.0_RKIND
rkSSPweights(4) = 0.0_RKIND
else
err = 1
call mpas_log_write('config_rk3_stages must 3 or 4', &
messageType=MPAS_LOG_ERR)
endif
else
err = 1
call mpas_log_write('config_time_integration = ' // trim(config_time_integration) &
Expand Down Expand Up @@ -299,8 +317,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
call mpas_timer_stop("advect thickness and tracers")
! If using SSP RK, then update thickness and tracers incrementally.
! For first RK stage, thickness and tracer updates above are sufficient
if ( rkStage > 1 ) then
! For first RK stage, thickness and tracer updates above are sufficient.
! Likewise, for the 4-stage SSP RK3 the last stage is just a forward euler update.
if ( (rkStage > 1) .and. (rkStage < 4) ) then
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d)
Expand Down Expand Up @@ -507,15 +526,16 @@ subroutine prepare_advection(domain, err)
real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity
real (kind=RKIND), pointer :: calvingCFLdt, faceMeltingCFLdt
integer, pointer :: processLimitingTimestep
integer, pointer :: processLimitingTimestep, config_rk_order, config_rk3_stages
integer, dimension(:), pointer :: edgeMask

logical, pointer :: config_print_thickness_advection_info
logical, pointer :: config_adaptive_timestep
logical, pointer :: config_adaptive_timestep_include_DCFL

character (len=StrKIND), pointer :: &
config_thickness_advection ! method for advecting thickness and tracers
config_thickness_advection, & ! method for advecting thickness and tracers
config_time_integration

integer :: &
allowableAdvecDtProcNumberHere, &
Expand Down Expand Up @@ -567,6 +587,9 @@ subroutine prepare_advection(domain, err)
call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep)
call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep_include_DCFL', config_adaptive_timestep_include_DCFL)
call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection)
call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration)
call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order)
call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages)

if (trim(config_thickness_advection) == 'none') then
if (config_adaptive_timestep) then
Expand Down Expand Up @@ -612,7 +635,6 @@ subroutine prepare_advection(domain, err)
err = ior(err, err_tmp)

allowableAdvecDtOnProc = min(allowableAdvecDtOnProc, allowableAdvecDt)

! Calculate diffusive CFL timestep, if needed
! This used to be only calculated if (config_adaptive_timestep_include_DCFL) but for simplicity,
! now it is always calculated. That allows assessment of the DCFL even when it is not being obeyed
Expand All @@ -637,6 +659,12 @@ subroutine prepare_advection(domain, err)
block => block % next
end do

! If using 4-stage SSPRK3, CFL number of 2 is theoretically allowed
if ( (trim(config_time_integration) == 'runge_kutta') .and. &
(config_rk_order == 3) .and. (config_rk3_stages == 4) ) then
allowableAdvecDtOnProc = allowableAdvecDtOnProc * 2.0_RKIND
allowableDiffDtOnProc = allowableDiffDtOnProc * 2.0_RKIND
endif

! Local advective CFL info
call mpas_set_timeInterval(allowableAdvecDtOnProcInterval, dt=allowableAdvecDtOnProc, ierr=err_tmp)
Expand Down

0 comments on commit 4f69cfe

Please sign in to comment.