Add embedded Runge-Kutta methods for automatic time-step adjustment#170
Open
trhille wants to merge 19 commits into
Open
Add embedded Runge-Kutta methods for automatic time-step adjustment#170trhille wants to merge 19 commits into
trhille wants to merge 19 commits into
Conversation
Also remove mask update from velocity solver and write to logs when updating and resetting masks.
Deallocate edgeMaskPrev and vertexMaskPrev where cellMaskPrev is deallocated.
Remove outdated comment about updating masks before velocity solve.
Add call to deallocate fluxAcrossGroundingLineOnCellsAccum
Control mask updates before velocity solve with a logical flag in the call to li_velocity_solve
Add RK pairs for embedded SSP RK time integration. These pairs were suggested by Co-pilot and are not necesarily aligned with pairs suggested in the literature.
Only allocate extra arrays used in embedded methods when embedded RK is actually being used. This will save memory overhead.
Synchronize embedded norm across ranks so we don't get a crash when different ranks have different norms, and thus each set a different dt.
Ignor damage and (nearly) ice-free regions for error norm calculation.
Make clock advancement consistent with new deltat when deltat is set by the embedded RK error norm. This means moving the clock update to a different point in the time step, which may cause changes relative to previous behavior.
Fix issue in set_timestep that was causing crash when using dt determined by embedded RK.
Floord RK-determined dt before advancing clock, which is necessary to avoid throwing an error with non-integer interval to next forcing time.
Use global rms error instead of max, which was causing dt to go to zero.
Apply face-melt speed after advection, so it uses the dt determined by the embedded RK method.
Use a SSP RK(4,2) for the embedded pair to RK(4,3), which prevents the time step from going to 0.
Use the optimized embedded pair from Fekete et al. (2022) for RK2, and the second-stage endpoint as a pair for RK(3,3). These work much better than the previous pairs; for a 4km Thwaites mesh with historical ISMIP6 forcing and 2nd/3rd order advection, they result in a small but stable time step of about 1.5 days. Meanwhile, the pair for RK(4,3) is allowing multi-month time steps, so it is unclear why these are so much more restrictive.
trhille
force-pushed
the
mali/embedded_ssp_rk_pairs
branch
from
e42b77a to
f581f11
Compare
There was a problem hiding this comment.
Pull request overview
This PR extends the MPAS-Albany Land Ice forward-mode time integration to support an embedded SSPRK error estimate with retry-based timestep control, while also adjusting when masks and grounded face-melt updates are applied to better preserve RK-stage consistency.
Changes:
- Add namelist options to enable/parameterize embedded SSPRK error estimation and retry-based timestep adjustment.
- Refactor the RK integrator to support step retries, deferred clock advance, and deferred grounded face-melt application until after step acceptance.
- Add an
updateMaskcontrol to velocity solves to avoid recalculating masks inside RK loops, and add additional mask-related logging.
Reviewed changes
Copilot reviewed 6 out of 6 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
| components/mpas-albany-landice/src/shared/mpas_li_mask.F | Adds a log line during mask updates. |
| components/mpas-albany-landice/src/Registry.xml | Introduces new namelist options for embedded SSPRK error control. |
| components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F | Adds updateMask flag to control mask updates before velocity solves. |
| components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F | Implements embedded SSPRK retry logic, timestep adaptation, and defers clock/melt application until acceptance. |
| components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F | Adds “compute vs apply” controls for grounded face melt to support deferred application. |
| components/mpas-albany-landice/src/mode_forward/mpas_li_core.F | Updates initial velocity solve call signature to pass updateMask. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| ! Calculate cellMask values=========================== | ||
| ! ==== | ||
|
|
||
| call mpas_log_write("Updating masks") |
Comment on lines
141
to
149
| 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 | ||
|
|
Author
There was a problem hiding this comment.
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.
Comment on lines
+1669
to
+1673
| if (.not. doApplyFrontAblation) then | ||
| block => block % next | ||
| cycle | ||
| endif | ||
|
|
Author
There was a problem hiding this comment.
The CFL calculation might be a legitimate concern, but we definitely do not want to calculate faceMeltingThickness yet.
Comment on lines
1674
to
+1690
| @@ -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) | |||
|
|
|||
| 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." |
Author
There was a problem hiding this comment.
Yes, it bothered me that co-pilot created an absolute tolerance field that is applied to multiple tracers with different units and magnitudes.
Remove declaration of unused variable Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
Weight error calculation in embedded RK method by cellArea
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This merge adds embedded pairs as error estimators for automatic, error-informed time-step adjustment when using Runge-Kutta time stepping. Each choice of Runge-Kutta scheme now comes with its own embedded pair, which is a lower-order scheme with the same number of stages. The difference between the two solutions is used as an estimate of the time discretization error. Relative and absolute error thresholds are set by the user in namelist, and the code accepts or rejects a time step based on the normalized error between the higher- and lower-order solutions. If a time step is rejected, the algorithm restarts the advection step with a smaller time step until the solutions converge, the specified allowable number of attempts is reached, or the time step goes below 1s.
Note that this also requires moving the application of face-melting and the clock advancement to be after advection. This is because the length of the time step is not known until the advection step has been accepted. There maybe be additional complications that result from this re-shuffle. The face-melt rate is still calculated at the same place in the time step, but it is not applied until we know the time-step length.
See Fekete et al. (2022) for an explanation of embedded methods. Note that we use their optimized weights for the pair to RK2, but we use different choices for RK(3,3) and RK(4,3), based on some trial and error. The optimized pairs from Fekete et al. (2022) for those schemes led to vanishing time-step size in the course of a few dozen time steps.
Also note that this contains the changes in #166.