Create a calving CFL limit and apply to adaptive timestepper

[matthewhoffman]

There is reason to think it is incorrect to allow calving to remove more than a whole grid cell in a given timestep - a CFL condition for calving. This merge calculates a calving CFL and adds it to the adaptive timestepper. It's important to note that the calving CFL is calculated from the previous timestep, because determining the calving CFL prior to starting the timestep is prohibitively difficult. That would require doing advection and other operations to generate the geometry that the calvingVelocity will be based on just to calculate the CFL limit, and then going back and actually progressing through the time step for real once the dt has been set. So this addition works under the principle that the calving CFL limit will not change a lot from timestep to timestep, thus a lagged CFL limit is a useful way to adapt the timestep. Testing shows this assumption works very well most, though not all, of the time.

The calving CFL is included in the adaptive timestepper with the new config_adaptive_timestep_include_calving option. The fraction of the calving CFL to use is controlled with the new option config_adaptive_timestep_calvingCFL_fraction.

A number of other related changes are made:

• The adaptive timestepper logic and log messages are adjusted somewhat. A significant change is that config_adaptive_timestep_CFL_fraction now only applied to the advective CFL. There is no longer a general fraction applied.
• The new variable processLimitingTimestep outputs a code indicating which process (advective CFL, diffusive CFL, calving) was the limiting factor in determining each timestep length.
• The new variable calvingCFLdt calculates the dt estimated for by the calving CFL. The variable dtCalvingCFLratio is the ratio of the actual dt used to that dt indicated by the lagged calving CFL calculation. If conditions don't change between timesteps, this should be the value chosen for config_adaptive_timestep_calvingCFL_fraction, but if calving conditions change significantly between timesteps, this variable can be used to assess how bad the choice of dt actually was.
• The 10% calving error now can be controlled with a namelist option: config_calving_error_threshold
  Note could add a check and a warning if the calving CFL on the actual time step exceeds the timestep in use.
• There was a longstanding problem in the adaptive timestepper of dying if the force interval required a timestep smaller than the minimum allowable timestep. This has been fixed.
• Log output has been modified to make the start of each timestep more obvious.

[matthewhoffman]

Results of full_integration suite:

Test Runtimes:
00:10 PASS landice_dome_2000m_sia_restart_test
00:05 PASS landice_dome_2000m_sia_decomposition_test
00:08 PASS landice_dome_variable_resolution_sia_restart_test
00:04 PASS landice_dome_variable_resolution_sia_decomposition_test
00:31 PASS landice_enthalpy_benchmark_A
00:13 PASS landice_eismint2_decomposition_test
00:13 PASS landice_eismint2_enthalpy_decomposition_test
00:14 PASS landice_eismint2_restart_test
00:14 PASS landice_eismint2_enthalpy_restart_test
00:18 FAIL landice_greenland_sia_restart_test
00:09 FAIL landice_greenland_sia_decomposition_test
00:09 PASS landice_hydro_radial_restart_test
00:06 PASS landice_hydro_radial_decomposition_test
00:17 PASS landice_dome_2000m_fo_decomposition_test
00:14 PASS landice_dome_2000m_fo_restart_test
00:11 PASS landice_dome_variable_resolution_fo_decomposition_test
00:13 PASS landice_dome_variable_resolution_fo_restart_test
00:14 PASS landice_circular_shelf_decomposition_test
00:41 FAIL landice_greenland_fo_decomposition_test
00:41 FAIL landice_greenland_fo_restart_test
00:32 PASS landice_thwaites_decomposition_test
00:44 PASS landice_thwaites_restart_test
00:51 FAIL landice_humboldt_mesh-3km_restart_test_velo-fo_calving-von_mises_stress_damage-threshold_faceMelting
Total runtime 07:14
FAIL: 5 tests failed, see above.

The greenland failures are because this PR applies the calving routine on the initial time, whereas before it was not applied until the first time step. For rate-based calving laws, applying calving on the initial time has no effect because dt=0, but for threshold-based calving laws, this change will apply calving on the initial time. The greenland cases all are currently using thickness-based calving, so applying calving on the initial time removes floating ice from a handful of grid cells in NE Greenland. As I mention in the commit messages, I think we actually prefer to have threshold-based calving applied on the initial time, so this change is desireable.

The landice_humboldt_mesh-3km_restart_test_velo-fo_calving-von_mises_stress_damage-threshold_faceMelting failure is due to a small change (order 1e-6 m/s) to calvingVelocity (but not calvingThickness). This PR includes a lot of minor adjustments to how calving is processed, so that level of change did not concern me.

Therefore, I would consider all these changes as blessed.

To actually test this PR, I would run the calving_dt_convergence suite in the accompanying COMPASS PR: MPAS-Dev/compass#386

@trhille , you may also want to run one of your Humboldt paper runs to ensure there are no major differences.

[trhille]

I see the same 5 tests failing baseline comparison.

Test Runtimes:
00:09 PASS landice_dome_2000m_sia_restart_test
00:05 PASS landice_dome_2000m_sia_decomposition_test
00:08 PASS landice_dome_variable_resolution_sia_restart_test
00:05 PASS landice_dome_variable_resolution_sia_decomposition_test
00:33 PASS landice_enthalpy_benchmark_A
00:12 PASS landice_eismint2_decomposition_test
00:12 PASS landice_eismint2_enthalpy_decomposition_test
00:14 PASS landice_eismint2_restart_test
00:14 PASS landice_eismint2_enthalpy_restart_test
00:18 FAIL landice_greenland_sia_restart_test
00:10 FAIL landice_greenland_sia_decomposition_test
00:09 PASS landice_hydro_radial_restart_test
00:06 PASS landice_hydro_radial_decomposition_test
00:16 PASS landice_dome_2000m_fo_decomposition_test
00:15 PASS landice_dome_2000m_fo_restart_test
00:11 PASS landice_dome_variable_resolution_fo_decomposition_test
00:13 PASS landice_dome_variable_resolution_fo_restart_test
00:14 PASS landice_circular_shelf_decomposition_test
00:41 FAIL landice_greenland_fo_decomposition_test
00:40 FAIL landice_greenland_fo_restart_test
00:32 PASS landice_thwaites_decomposition_test
00:43 PASS landice_thwaites_restart_test
00:51 FAIL landice_humboldt_mesh-3km_restart_test_velo-fo_calving-von_mises_stress_damage-threshold_faceMelting
Total runtime 07:13
FAIL: 5 tests failed, see above.

[trhille]

I also tested this on the Humboldt_1to10km_r04_20210615 mesh using the most aggressive mass loss scenario from the perturbed parameter ensemble in my paper (HadGEM2 climate, 170 kPa von Mises threshold stress, q = 1/7). Happily, these changes result in only minor differences from the original simulation, which used an executable compiled using https://github.com/MALI-Dev/E3SM/tree/d0c286778a6a035c600af891a268ed9f843aac52 and https://github.com/sandialabs/Albany/tree/ecbc0182015fda8d6a2bb29eaa3a47c260f40b13.
In the results below, blue is the run from the original ensemble, green is compiled on the current head of MALI-Dev/develop (commit 6c06533, with commit d4f7c0b cherry-picked to enable compiling on Cori), and orange was compiled using this branch (commit 258eadc, also with commit d4f7c0b cherry-picked):

The green curve used config_distribute_unablatedVolumeDynCell = .true., but ran into the >10% ablation error, and I did not bother reducing the timestep and restarting.

The orange curve (this branch) used

config_adaptive_timestep_calvingCFL_fraction = 1.0
config_adaptive_timestep_include_calving = .true.
config_distribute_unablatedVolumeDynCell = .false.

The ratio of timestep taken to CFL calculated at the end of the timestep (i.e., the effect of the one-timestep lag in CFL condition calculation) looks pretty good here:

[trhille]

I ran the calving_dt_convergence COMPASS test suite (MPAS-Dev/compass#386), and all tests passed successfully:

Test Runtimes:
02:42 ESC[92mPASSESC[0m landice_calving_dt_convergence_mismip+.specified_calving_velocity.none
02:12 ESC[92mPASSESC[0m landice_calving_dt_convergence_mismip+.von_Mises_stress.none
02:01 ESC[92mPASSESC[0m landice_calving_dt_convergence_humboldt.specified_calving_velocity.none
01:48 ESC[92mPASSESC[0m landice_calving_dt_convergence_humboldt.von_Mises_stress.none
03:03 ESC[92mPASSESC[0m landice_calving_dt_convergence_thwaites.specified_calving_velocity.none
03:28 ESC[92mPASSESC[0m landice_calving_dt_convergence_thwaites.von_Mises_stress.none
73:56 ESC[92mPASSESC[0m landice_calving_dt_convergence_mismip+.von_Mises_stress.FO
08:40 ESC[92mPASSESC[0m landice_calving_dt_convergence_humboldt.von_Mises_stress.FO
21:36 ESC[92mPASSESC[0m landice_calving_dt_convergence_thwaites.von_Mises_stress.FO
Total runtime 119:29
PASS: All passed successfully!

See MPAS-Dev/compass#386 for test results.

[trhille]

@matthewhoffman, this is great! I had a few minor suggestions here and there. One concern that could complicate things is that I think we should include face-melting in the CFL condition calculation, although it is unlikely to exceed calving rates.

[matthewhoffman]

@trhille , thanks for the review and the informative testing on this. I updated the PR. Let me know if you think we should add the faceMelting CFL limit or not. I'm inclined to not bother, but I could be convinced.

[trhille]

@matthewhoffman, thanks for these changes. I don't think any of them will make a difference, but I'll run full_integration again to be sure.

@trhille , thanks for the review and the informative testing on this. I updated the PR. Let me know if you think we should add the faceMelting CFL limit or not. I'm inclined to not bother, but I could be convinced.

Let's leave it as-is for now. It's probably true that in the vast majority of cases, face melting will not be the dominant velocity, and for our short-term Antarctic use cases that won't matter. We can add that in another PR if we determine it is necessary.

[trhille]

Okay, I'm seeing the same behavior as before when I run full_integration, so this is good to go!