Update weights for ocean barotropic subcycling in split-explicit solver

[mark-petersen]

These improvements were originally written by @dengwirda, as described in E3SM-Ocean-Discussion#48. Please refer to that PR for additional test results and discussion. Also see test results run in combination with AB2 at #5989 (comment).

This PR updates the weights config_btr_gam1_velWt1, config_btr_gam2_SSHWt1 in the MPAS-Ocean barotropic solver based on recent analysis of this scheme. This update applies to "split explicit" time stepping schemes, i.e. config_time_integrator = 'split_explicit' and the new config_time_integrator = 'split_explicit_ab2' in #5989. The new weights allow for a barotropic time step config_btr_dt to be 25% longer in the tests of EC30to60, thus speeding up the barotropic subcycling. We expect that values of config_btr_dt can be increased by 25% for all meshes.

[NML]
[non-BFB]

[mark-petersen]

@dengwirda, @jeremy-lilly, @gcapodag: Many thanks for your work on time-stepping methods this past year, and this tuning in particular. Feel free to comment here as time permits.

[mark-petersen]

Adding previous comment for completeness (and the great images!!!):

This PR updates the default gamma weights (config_btr_gam1_velWt1, config_btr_gam2_SSHWt1) in the barotropic subcycling phase of the split-explicit solver. Based on work with @jeremy-lilly, @gcapodag and @mark-petersen, a stability analysis of the linearised external mode problem was used to estimate the coefficients that improve the scheme's CFL-limit, while also reducing numerical dissipation. Overall, the new coefficients $\gamma_{1,2} = \frac{8}{15}$ allow larger time-steps to be used; approximately $\frac{4}{3}$ the length of the maximum stable step-size used with the current $\gamma_{1} = \frac{1}{2}$, $\gamma_{2} = 1$ pair. This allows, e.g. config_btr_dt to be increased from 60sec to 75sec for QU30km and EC60to30km meshes and overall runtime to be reduced.

Following Shchepetkin et al., the behaviour of a time-stepping scheme can be analysed using a dissipation-dispersion diagram, in which a 'perfect' integration scheme would reproduce the unit-circle exactly, stability requires the scheme to lie inside the circle, dissipative errors are represented by a reduced curve radius and dispersive errors by angular offset.

Old weights:

New weights:

The current weights achieve $\mathrm{CFL_{max}} = 1.41$ and the new weights $\mathrm{CFL_{max}} = 1.88$, leading to the $\frac{4}{3}$ time-step extension.

Though there is a third weight config_btr_gam3_velWt2 available, this is typically not used, with config_btr_solve_SSH2 = .false. Initial experiments did not suggest including this step improved either stability or accuracy (the opposite in fact), though additional investigation would be required on this front.

I've confirmed the new weights allow config_btr_dt = 00:00:75 for QU30km and EC60to30km config.'s (vs the current 60 sec time-step) and that 3yr standalone spin-ups lead to eyeball-norm consistent (but non-BFB) output. Additional testing would be required to validate performance in coupled E3SM runs.

• Shchepetkin, A.F. and McWilliams, J.C., 2009. Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations. In Handbook of Numerical Analysis (Vol. 14, pp. 121-183).

[mark-petersen]

To take full advantage of this PR, we should lengthen config_btr_dt where we can for other meshes. @hyungyukang could you please extend your 18to6 run a bit (10 days?) with the new btr weights, and see if you can use config_btr_dt of 15s (or even 16?) rather than 12s?

[hyungyukang]

@mark-petersen , Sure, I'm going to run the same 18to6 run I've done with updated btr weights and see if we can use larger BTR time step size. I will report here once I finish the tests.

[mark-petersen]

@jonbob, @darincomeau and @xylar, I suspect these three could all be extended from 15s to 19 or 20s. I think it is too much to test all this right now for this PR, but as we test and add E3 version grids, we should test with the value of config_btr_dt that is 25% longer.

[mark-petersen]

I ran one year of this PR plus #5989 with gnu on chrysalis in MPAS-O standalone, using compass ec30to60 and ecwisc30to60 suites, but with cached mesh and init steps to save time. This then takes the model through the spin up process, then I extended the simulation step to 360 days. The results for AB2 look very good. It is stable with the same spin-up procedure, and min and max on temperature and salinity remain constant.

[dengwirda]

@mark-petersen and all, as you've said, these changes should allow config_btr_dt to be lengthened in general, and I previously had success across a range of meshes + standalone spin-ups wrt. similar % increases --- would expect similar impact irrespective of resolution, etc. Of course, tuning for max config_btr_dt will be a per-mesh experience, so those tests would have to be done for each config. individually...

In addition to the longer config_btr_dt, these new weights should also lessen (spurious) dissipation of gravity waves. This was observable in a channel config., though I don't expect this behaviour would feedback significantly onto the mean climate state. @jonbob previously confirmed the expected small diffs in a 30-to-60km B-case (G-case?).

Overall, hopefully this is simply some 'free' MPAS-O performance improvement, via optimisation of the time integrator CFL + truncation behaviour.

[rljacob]

@hyungyukang please review.

[hyungyukang]

@rljacob and @mark-petersen , sorry for my late review.
I ran some cases but forgot to comment here.

I ran several real-world test cases on RRS18to6v3 mesh using the standalone MAPS-O with updated weights and tried to increase btr dt from 12 s (default). Both AB2 and forward-backward schemes were stable for 18 s of btr dt, where AB2 uses 75% of the default del 4 coefficient (config_mom_del4=2.4e09) and config_mom_del4_div_factor=4.0. Below figure shows the average and maximum kinetic energy. Both schemes were blown-up when btr dt is larger than 18 s. This result can be different for other high-resolution simulations, but at least 15 seconds seems stable on the RRS18to6 mesh.

[hyungyukang]

I'm approving this PR based on visual inspection and testing from @dengwirda , @mark-petersen , and @hyungyukang (also in #5989).

[jonbob]

passes sanity testing with expected DIFFs -- merged to next

[jonbob]

merged to master and expected DIFFs blessed