Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'trhille/implement_rk_time_integration' into MALI-Dev/de…
…velop 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/implement_rk_time_integration: (36 commits) Move subglacial hydro and bed topography updates back to their former positions Always update velocity after calving if it was not updated before calving. Make config_update_velocity_before_calving = .false. the default Add logic to control whether velocity is solved before and/or after calving Average fluxAcrossGroundingLineOnCells in time integration module Revert "Include calving in Runge-Kutta loop" Revert "Calculate layerThickness after calving in RK loop" Calculate layerThickness after calving in RK loop Include calving in Runge-Kutta loop Remove support for config_rk_order = 1 Restore calving front within RK loop Initialize RK weights to large negative values Simple clean-up after code review Update temperature and waterFrac in RK loop, but not enthalpy Update thickness and tracer halos before velocity solve Update stresses and strain rates after calving Reset layerThickness when using config_restore_thickness_after_advection Calculate groundedToFloatingThickness and grounding line flux after RK Fix bug in groundedSfcMassBalApplied after Runge-Kutta loop Some more cleanup ...
- Loading branch information
