Implement Runge-Kutta time integration

[trhille]

This merge adds Runge-Kutta time integration to MALI, which is likely necessary when using higher-order advection. Currently, Runge-Kutta integration is used for damage evolution and thickness and tracer advection. Surface and basal mass balances are applied within the Runge-Kutta loop and budgets are calculated at the end of the loop. Front ablation (i.e., calving and face-melting) and bed topography updates are applied via forward Euler outside the Runge-Kutta loop. This could conceivably be updated for full consistency.

The Runge-Kutta formulations used here are the two-stage second-order and and three-stage third-order strong stability preserving schemes of Shu and Osher (1988), as well as the four-stage third-order strong stability preserving scheme of Spiteri and Ruuth (2002). See equations 2.47–2.49 in Durran (2010) for an overview.

There is one velocity solve per Runge-Kutta stage and an optional solve before calving (more on that below), meaning that the four-stage RK3 scheme is ~25–33% more expensive than the three-stage RK3 scheme for a given time step length. However, the four-stage scheme theoretically allows for an effective CFL fraction of 2.0, while the three-stage scheme is limited to CFL fraction of 1.0. Testing indicates stable results using the four-stage scheme for an effective CFL fraction of 1.8 for a simple grounded slab advecting at uniform speed. Note that when using the four-stage scheme in MALI, the maximum allowable time step is updated, so config_adaptive_timestep_CFL_fraction should still be between 0.0 and 1.0.

There is now an optional extra velocity solve between advection and calving in addition to the final velocity solve at the end of the time step, for any choice of time integration scheme. The extra velocity solution prior to calving ensures a self-consistent set of inputs to calving routines; i.e., the velocity, stress, and strain rate fields will be consistent with the geometry used for calving, which is not the case when velocity is solved only after calving. This makes calving more accurate, but is obviously significantly more expensive. The extra velocity solution is enabled using config_update_velocity_before_calving = .true.. It is disabled by default to allow regression tests to pass baseline comparison. Note that the extra velocity solution is not necessary for some calving laws, such as damage threshold calving, so the user should carefully consider this option when setting up simulations. There is currently no check in the code to ensure the combination of calving and extra velocity solve settings are sensible. Also of note is that li_restore_calving_front (if enabled) is called before the extra velocity solution to prevent solving velocities on an extended geometry, while the remaining calving routines are called after the velocity solve. If the volume of dynamic ice is not changed by calving, then the final velocity solution is automatically skipped.

References:
Durran, Dale R. Numerical Methods for Fluid Dynamics. Vol. 32. Texts in Applied Mathematics. New York, NY: Springer New York, 2010. https://doi.org/10.1007/978-1-4419-6412-0.

Shu CW, Osher S (1988) Efficient implementation of essentially nonoscillatory shock-capturing schemes. J Comp Phys 77:439–471

Spiteri RJ, Ruuth SJ (2002) A new class of optimal high-order strong-stability-preserving time discretization methods. SIAM J Numer Anal 40:469–491

[trhille]

@matthewhoffman, I expect that this will require significant refactoring for flexibility, but I think having the in-progress PR open now will be useful for discussion and testing. Note that only the last six commits here are really relevant, as the others are due to merging my FCT branch into this one. I can clean up that history once the FCT branch is merged into MALI-Dev/develop.

[trhille]

One outstanding issue: what do we do about calculating fluxAcrossGroundingLine?
https://github.com/trhille/E3SM/blob/implement_rk_time_integration/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F#L630-L659

It seems like this should be moved to occur after the RK loop. This is also a first-order treatment. Is it possible to update it to higher order when using fct advection?

[trhille]

Testing

I ran the landice/mesh_convergence/horizontal_advection case (MPAS-Dev/compass#690) for 8, 4, 2, and 1km resolutions with the following namelist settings:

config_thickness_advection = 'fct'
config_tracer_advection = 'fct'
config_horiz_tracer_adv_order = 3
config_advection_coef_3rd_order = 1.0
config_time_integration = 'runge_kutta'
config_rk_order = 3
config_rk3_stages = 3

We get close-to-third-order convergence for third order advection with three-stage RK3:

Using the same settings but with config_rk3_stages = 4 for four-stage RK3, and setting time step to dt*2 in forward.py:

Using second-order advection and RK2:

config_thickness_advection = 'fct'
config_tracer_advection = 'fct'
config_horiz_tracer_adv_order = 2
config_advection_coef_3rd_order = 1.0
config_time_integration = 'runge_kutta'
config_rk_order = 2

[trhille]

To check that SMB budget conserves mass, I ran the grounded slab test case with pure third-order advection and a uniform SMB of 1 m/yr.
Plots below are difference between grounded mass change and SMB for each RK case.
Three-stage RK3:

Four-stage RK3:

RK2:

These were calculated changing one line in plot_mass_balance.py

# get fields and plot stuff
volGround = f.variables['groundedIceVolume'][:]
volGround = volGround * rhoi / 1.0e12
dvolGround = np.zeros(yr.shape)
dvolGround[1:] = (volGround[1:]-volGround[:-1]) / dyr[1:]
#ax.plot(yr, dvolGround, 'k', linewidth=3, label="grounded mass change")

SMBg = f.variables['totalGroundedSfcMassBal'][:] / 1.0e12
ax.plot(yr, dvolGround - SMBg, label="mass change minus grounded SMB")

[trhille]

Tested mass conservation using both surface and grounded basal mass balance. SMB treated same as above. For BMB, I set basalHeatFlux = 1.0 everywhere and added:

    config_thermal_solver = 'temperature'
    config_thermal_calculate_bmb = .true.

Three-stage RK3:

Four-stage RK3:

RK2:

[trhille]

Results of slotted cylinder test case using config_horiz_tracer_adv_order = 3; config_advection_coef_3rd_order = 0.5:
CFL fraction 0.99: https://drive.google.com/file/d/1TdydxWOpMwf76sOhQrbiTp_rBluyuTlo/view?usp=sharing
CFL fraction 0.25: https://drive.google.com/file/d/1kLvIexVLpmyvIeD5ye_7ngktQlrJN2yV/view?usp=sharing

[trhille]

One outstanding issue: what do we do about calculating fluxAcrossGroundingLine? https://github.com/trhille/E3SM/blob/implement_rk_time_integration/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F#L630-L659

It seems like this should be moved to occur after the RK loop. This is also a first-order treatment. Is it possible to update it to higher order when using fct advection?

Updating this to higher order in #94. Once that is merged I will rebase this branch and try to calculate the final fluxAcrossGroundingLine for the time step as a weighted average, in the same way as the surface and basal mass balance terms. I think that should close the budget, but we'll have to try it to find out.

[trhille]

Testing the calculation of groundedToFloatingThickness and fluxAcrossGroundingLine after the Runge-Kutta loop (introduced in 3ec981a). Currently this will only be correct for fo thickness advection, as higher order grounding line flux (#94) is awaiting review. I compared a run on this branch using both the 4-stage RK3 and forward Euler integration with a run compiled at commit 6d0d572 using forward Euler. Neither budget completely closes, likely due to the same issues we've been having with budgets for a long time (see PR #38 and Issue #23).

Forward Euler (6d0d572):

Forward Euler (3ec981a):

RK3 (3ec981a):

fractional error is calculated as (SMBg+BMBg-GLflux-GLMigrationflux - dvolGround)/dvolGround

[trhille]

Here are full_integration results as of 4005761:

Forward Euler: all tests pass execution and validation. I have not yet run a baseline comparison.
RK2: All tests pass execution, three fail validation:

• landice_hydro_radial_restart_test. waterThickness and waterPressure L_inf norms of 1.5e-15 and 4.7e-10, respectively
• landice_dome_2000m_fo_decomposition_test. thickness L2 norm of 1.7e-11 and normalVelocity L_inf norm of 1.1e-18
• landice_dome_variable_resolution_fo_decomposition_test. normalVelocity L1 and L2 norms of 8.8e-17 l2: 1.68e-18

3-stage RK3: All tests pass execution, one fails validation:

• landice_hydro_radial_restart_test. waterThickness and waterPressure L_inf norms of 1.7e-15 and 4.7e-10, respectively

4-stage RK3: All tests pass execution, one fails validation:

• landice_hydro_radial_restart_test. waterThickness and waterPressure L_inf norms of 1.3e-15 and 4.7e-10, respectively

[trhille]

A few data points on the hydro_radial_restart test that may or may not be helpful:

1. Moving SGH solve after last RK stage (but before velocity solve) did not help.
2. Updating halos on cellMask, edgeMask, and vertexMask before calling SGH solve did not help.
3. Changing config_thermal_calculate_bmb = .true. to .false. did not help.
4. Using groundedBasalMassBalApplied instead of groundedBasalMassBal within the hydro module did not help.
5. The test fails validation on both Perlmutter and Chicoma.
6. Increasing number of halos from 3 to 4 did not help.

[matthewhoffman]

The small differences for the landice_dome_2000m_fo_decomposition_test and landice_dome_variable_resolution_fo_decomposition_test are no concern, and we can just loosen the tolerances for those tests by an order of magnitude. Thanks for digging in to the hydro failures. Those are perplexing given the hydro solver is not directly affected by the time-stepping method, and nothing outside the hydro model is evolving in these tests. I'll look at the code to see if I have any ideas about what could be going on there.

[matthewhoffman]

@trhille , the level of documentation you've made in this PR is really great. A few items for you:

1. What are you currently thinking about the treatment of restore_calving in the RK loop? Having velocity solved on an expanded ice extent in the case when we impose a fixed ice extent seems potentially problem-causing. That said, I'd also be happy to leave that for a separate PR as the RK and FCT methods are not yet used in any production configurations.
2. In the description, mention that the FE treatment is unified with the RK treatment, and therefore the order of operations are changed for FE as well.

[matthewhoffman]

@trhille , as monumental as this PR is, the implementation seems fairly straightforward, and I only have minor comments/requests. The big work is probably in evaluating and blessing the changes to all the compass test cases. Have you looked at baseline comparisons for Forward Euler yet? I wonder how gruesome it is. Given this PR is primarily bringing in new functionality that is not currently in production runs, and you've demonstrated that functionality is working as expected, at least to our current expectations, if the changes from baseline tests with FE make sense and are palatable, I don't see any obstructions to getting this merged in the next week or so.

Regarding the hydro_radial restart weirdness, I'm not entirely opposed to merging this without resolving that, but maybe after having looked this over carefully I'll have other ideas about what to look into.

[trhille]

Thanks for the review, @matthewhoffman! I'm glad it was more straightforward to review than we probably both anticipated.

Regarding the hydro_radial restart weirdness, I'm not entirely opposed to merging this without resolving that, but maybe after having looked this over carefully I'll have other ideas about what to look into.

It's definitely tempting to merge without fixing the hydro_radial restart issue, but it also makes me nervous because then we have a baseline that we know to be incorrect for that test.

[matthewhoffman]

@trhille , do you have any results on relative computational cost between these schemes? (e.g. compass timers) I don't want a whole performance study for this PR, but a rough sense of the impact would be good to know. (Though that may change depending on how we handle the 'extra' calving velocity solve issue.)

[trhille]

Compass runtimes on Chicoma:
Forward Euler: 16:26
RK2: 18:40
3-stage RK3: 20:25
4-stage RK3: 21:31

Note that landice_hydro_radial_decomposition_test now takes about 6.5 minutes on Chicoma, for some reason, when it should only take about 6 seconds.

[trhille]

1. What are you currently thinking about the treatment of restore_calving in the RK loop?  Having velocity solved on an expanded ice extent in the case when we impose a fixed ice extent seems potentially problem-causing.  That said, I'd also be happy to leave that for a separate PR as the RK and FCT methods are not yet used in any production configurations.

@matthewhoffman, sorry I somehow missed this question. This is a good point. We could include a call to li_restore_calving_front within the RK loop before the velocity solve, but deal with any other calving outside of the RK loop. Or we could accept that for purposes of RK integration the calving front may be allowed to advance beyond the initial extent with the knowledge that it will be dealt with before the end of the time step. If we decide to add yet another velocity solve after calving, that might be the way to go. But as you say, this might be better as the subject of another PR.

Update: I've move the call to li_restore_calving_front to within the RK loop prior to the velocity solve in bd3d61e.

[trhille]

Testing mesh convergence in a simplified test case with calving at e2a2f0b. I used the circular_shelf COMPASS test to create four shelves with 40 km radius at 10km, 5km, 2.5km, and 1.25km resolution. I gave them very small, radially symmetric velocities:
ncap2 -s uReconstructX = xCell / (1.0e5* 3.154e7); uReconstructY = yCell / (1.0e5 * 3.154e7) landice_grid.nc

I then ran these forward for 10 years on one processor using the following relevant namelist settings, designed to result in a shelf of 30 km radius at the end of the run:

    config_velocity_solver = 'none'

    config_thickness_advection = 'fct'
    config_tracer_advection = 'none'
    config_horiz_tracer_adv_order = 3
    config_advection_coef_3rd_order = 1.0

    config_calving = 'specified_calving_velocity'
    config_calving_specified_source = 'const'
    config_calving_velocity_const = 3.17415e-5 ! 1km/yr

    config_rk_order = 3
    config_rk3_stages = 3
    config_adaptive_timestep = .true.
    config_min_adaptive_timestep = 0.0
    config_max_adaptive_timestep = 3.15e9
    config_adaptive_timestep_CFL_fraction = 0.25
    config_adaptive_timestep_calvingCFL_fraction = 0.25

I then created the 30km-radius shelves at the same resolutions and used these as the "true" solution for each run. I adapted the analysis script from the ocean/planar_convergence/horizontal_advection COMPASS test to determine convergence order.
I used both area and volume error metrics, which give slightly different convergence orders:

[trhille]

Regression Testing

full_integration results at e2a2f0b, using different config_time_integration options with the default COMPASS advection settings:
Forward Euler: all tests pass validation and baseline comparison. Run time on Chicoma = 15:56
NB: All run times include hydro_radial/decomposition_test taking over 6 minutes on Chicoma when it should take 6 seconds.
RK2: Run time on Chicoma = 18:16
Tests failing validation:

• landice/hydro_radial/restart_test
  • waterThickness: l1: 1.69114722226027e-13 l2: 8.09609254772241e-15 linf: 2.10942374678780e-15
  • waterPressure l1: 9.88184183370322e-08 l2: 5.73213291028481e-09 linf: 4.65661287307739e-10
  • NB: these are could be caused by a roundoff error in thickness, which is evident in the humboldt hydro case, below.
• landice/humboldt/mesh-3km_restart_test/velo-none_calving-none_subglacialhydro
  • strangely, there are thickness errors on the order of 1e-13m between the full and restart runs. This is probably causing all the other errors in waterThickness (O(1e-14),hydropotential (O(1e-9), etc.
• landice/dome/2000m/fo_decomposition_test
  • thickness: l1: 2.37392664553804e-10 l2: 1.71442952818835e-11
  • normalVelocity: linf: 1.14052986347792e-18
• landice/dome/variable_resolution/fo_decomposition_test
  • normalVelocity l1: 8.76103653375565e-17 l2: 1.62421391384178e-18

3-stage RK3: Run time on Chicoma = 20:00
Tests failing validation:

• landice/hydro_radial/restart_test
• landice/humboldt/mesh-3km_restart_test/velo-none_calving-none_subglacialhydro

4-stage RK3: Run time on Chicoma = 21:07
Tests failing validation:

• landice/hydro_radial/restart_test
• landice/humboldt/mesh-3km_restart_test/velo-none_calving-none_subglacialhydro

All tests using Runge-Kutta time integration fail baseline comparison except for landice_circular_shelf_decomposition_test, which has no time evolution.

This includes failure of baseline comparison for the landice/enthalpy_benchmark/A test, which is somewhat confusing because there is no geometry evolution. However, the errors are small enough that they may stem from round-off error in the reconstruction of the forcing and/or temperature fields in the RK loop.
Here are the error norms from the last time step of the 3-stage RK3 compared the the forward-Euler baseline:
Temperature: 149: l1: 5.80939740757458e-10 l2: 4.31429191711591e-11 linf: 3.58113538823090e-12
basalWaterThickness: 79: l1: 3.68393475058681e-07 l2: 1.84196737529341e-07 linf: 9.20983687646704e-08
groundedBasalMassBal: 79: l1: 6.14892980147419e-20 l2: 3.07446490073709e-20 linf: 1.53723245036855e-20

I also tested landice/humboldt/mesh-3km_restart_test/velo-none_calving-none_subglacialhydro using RK2 with config_SGH = .false.. This also failed validation, with the same thickness errors as before: 90: l1: 2.75890560397229e-11 l2: 2.08855832908240e-12 linf: 4.54747350886464e-13. Plotting up the differences between the first and last output times shows roundoff-level errors in both full_run and restart_run, but they are not the same. An equivalent run using forward Euler time stepping showed all zeros. So apparently there are two issues when using Runge-Kutta time stepping with advection turned off: round-off error in thickness and a non-BFB restart.

Upon further testing, it turns out that I can get the hydro restart tests to pass validation by adding a call to li_calculate_layerThickness again after layerThickness is summed into thickness, here: https://github.com/trhille/E3SM/blob/implement_rk_time_integration/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F#L387
However, that does not get rid of the roundoff error issue. It simply makes the roundoff errors consistent between the full run and the restart tun. Given that the case of no advection with higher-order time integration is basically useless, we have decided not to address this issue in this PR.

[trhille]

Mesh Convergence Testing

Testing mesh convergence using Matt's COMPASS branch at matthewhoffman/compass@080d1f0, using 1, 2 ,4, and 8km meshes. I used the following relevant namelist settings:

    config_thickness_advection = 'fct'
    config_tracer_advection = 'fct'
    config_horiz_tracer_adv_order = 3
    config_advection_coef_3rd_order = 1.0

    config_time_integration = 'runge_kutta'
    config_rk_order = 3
    config_rk3_stages = 3

horizontal_advection tests give close to third-order convergence:


For comparison, FO advection with forward Euler time integration yields less-than-first-order convergence:

Meanwhile, the Halfar dome test gives first order convergence. It's unclear if that is because of the SIA velocity, the first-order treatment at the ice margins, or something else. Running for 1000 years instead of 200 years did not improve convergence order.

[trhille]

Even more testing

I ran the fct_integration suite (trhille/compass@c884eaf) for the following cases: 2nd-order advection with RK2, 3rd-order FCT (beta = 0.5) with 3-stage RK3, and 3rd-order FCT (beta = 0.5) with 4-stage RK3. All tests passed validation, except for the 2nd-order landice_dome_variable_resolution_fo_fct_decomposition_test, which had the following error norms:
thickness: 2: l1: 1.54766681241536e-11 l2: 1.69923778717874e-12 linf: 4.54747350886464e-13
normalVelocity: 1: l1: 8.22396491678941e-17 l2: 1.52703232416953e-18 linf: 9.42324153820409e-20
which are small enough to be ignored.

[trhille]

Bonus convergence testing

2nd order advection, RK2:

3rd order advection, 4-stage RK3:

3rd order FCT, beta = 0.75, 3-stage RK3:

3rd order FCT, beta = 0.5, 3-stage RK3:

3rd order FCT, beta = 0.25, 3-stage RK3:

[matthewhoffman]

@trhille , this is amazing work! And the testing... I can't find anything to nitpick about and will merge this now. As we discussed, there are some subtleties with the RK treatment we don't fully understand, but all the meaningful configurations we want to support look great and the existing Forward Euler behavior is bit-for-bit. Thanks for documenting all the testing results.