Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add and revise step_MOM_dyn_split_RK2b #1

Conversation

theresa-morrison
Copy link

Directly copied from Bob's PR #534
"This pull request consists of two commits that add the new module MOM_dynamics_split_RK2b and then revise it to give a new split time-stepping scheme. This new scheme is enabled by setting the new runtime parameter SPLIT_RK2B = True.

With this new scheme, the velocity that is passed back and forth to the driver is averaged over the phases of the barotropic stepping, rather than being an instantaneous velocity. The instantaneous velocities can be regenerated from the phase-averaged velocities using a stored barotopic velocity increment. It is hoped that this new scheme will be slightly more robust than before, although this will be more expensive than the other (previous) split time stepping scheme, with two calls to horizontal_viscosity per baroclinic dynamics time step (instead of one) and four calls to continuity per step (instead of 3). Some of this added cost could be recovered by unwinding the sub-calls within continuity to reuse the PPM thickness reconstructions in one of the two directions.

This new scheme has been tested on all of the examples in the MOM6-examples test suite, where it gives very similar answers. It has also been confirmed that the new scheme reproduces across processor layout and restart files.

When this new scheme is used there are substantial differences in the contents of the restart files, with 7 fewer 3-d fields, but two additional 2-D velocity increments. Moreover, when this scheme is used in ALE mode, there is only a single pair of velocity components that need to be remapped.

Answers will change with SPLIT_RK2B = True. There is a new publicly visible module with which change and several new public interfaces. In addition, there are changes to the contents of the MOM_parameter_doc.all files even when the new scheme is not being exercised.

The commits in this PR include:

9feaa3b *+Revise algorithm in step_MOM_dyn_split_RK2b
c075fb8 +Add MOM_dynamics_split_RK2b enabled by SPLIT_RK2B"

We would like to test these changes in our regional domains ASAP, so are short cutting pulling from dev/gfdl only.

Hallberg-NOAA and others added 4 commits December 14, 2023 15:56
  Add the new module MOM_dynamics_split_RK2b and calls to the routines in this
module to MOM.F90.  These calls are exercised when the run time parameter
SPLIT_RK2B is true.  For now, all answers are bitwise identical when this new
code is being exercised, but there is a new module with multiple public
interfaces and a new entry in many MOM_parameter_doc files.
  Revised step_MOM_dyn_split_RK2b to use the time-filtered velocities as
arguments and reconstruct the instantaneous velocities with the proper phase
from the barotropic solver from the saved increments du_av_inst and dv_av_inst.
As a part of this change, the continuity solver is used to find the thickness
fluxes used in the predictor step Coriolis terms, and the horizontal viscous
accelerations are calculated for both the predictor and corrector steps, thereby
avoiding the need to vertically remap any 3-d fields apart from a single pair of
velocity components.

  The run-time option STORE_CORIOLIS_ACCEL is no longer used when SPLIT_RK2B =
True.  FPMIX was also disabled in this mode due to unresolved dimensional
inconsistency errors in that code.

  Additionally, the time-filtered velocities are now used instead of the
instantaneous velocities in the second call to radiation_open_bdry_conds(),
which should improve the performance of several of the Orlanski-type open
boundary conditions.

  The 2-d fields du_av_inst and dv_av_inst were added to the restart file, while
the 3-d fields u2, v2, diffu and diffv and either CAu and CAv or uh, vh and h2
are all removed from the restart files.  Remap_dyn_split_RK2b_aux_vars() now has
nothing to do, and it should probably be eliminated altogether in a subsequent
commit.

  All answers are changed when SPLIT_RK2B = True, and there are changes to the
contents of the restart files.
  This commit includes further revisions to MOM_dynamics_split_RK2b that avoid
some unnecessary calculations and group some of the halo updates into fewer
group passes.  All answers are bitwise identical.
alex-huth and others added 5 commits December 20, 2023 15:39
…d/flipped input and grids, and under dimensional rescaling to the +/-140 power. Bitwise identical results were achieved under rotation/rescaling using standard 2D test cases such as MISMIP+. Major changes include the following:

- Specified order of operations throughout, and modified the implementation of FEM integration,  h-point-to-node operations, etc for consistency under rotation (see MOM_ice_shelf_dynamics.F90 subroutines calc_shelf_driving_stress, CG_action, CG_action_subgrid_basal, matrix_diagonal, shelf_advance_front, and interpolate_h_to_b).
- Added an ice-shelf version of 'first_direction' and 'alterate_first_direction', allowing users to control the order in which the x- and y-direction ice-shelf advection calls are made. Required for consistency under rotation.
- Dimensional rescaling fixes throughout MOM_ice_shelf_initialize.F90, and within MOM_ice_shelf_dynamics.F90 (see subroutines initialize_ice_shelf_dyn, ice_shelf_solve_outer, and calc_shelf_taub)
- Rotation/rescaling testing revealed two additional bugs that were subsequently fixed:
(1) An index error in the conjugate gradient algorithm (subroutine CG_action)
(2) Discretization error for the east/north fluxes in subroutines ice_shelf_advect_thickness_x and ice_shelf_advect_thickness_y.

The commit also includes several simple changes for code readability and computational efficiency:

- Saved Phi, PhiC, and Phisub in the ice shelf dynamics control structure at the start of each run, rather than reallocating and recalculating these static fields each time they are used. Reshaped the Phisub array for computational efficiency over loops.
- Modified the sub-cell grounding scheme so that it is only called for partially-grounded cells that contain the grounding line (i.e. where float_cond(i,j)==1). Added float_cond to the ice dynamics structure for output.
- Simplified the option to calculate ice viscosity at 4 quadrature points per cell rather than at cell centers, by adding parameter NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS and reshaping CS%ice_visc accordingly.
- Style changes and removal of unused code (e.g. removed subroutine apply_boundary_values and field thickness_bdry_val)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants