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

Fix biharmonic leith #268

Merged
merged 9 commits into from
Feb 10, 2024
Merged

Conversation

iangrooms
Copy link

This PR fixes some bugs in Leith+E, changes the code so that it preserves rotational symmetry, re-orders some computations to improve efficiency, and introduces some new features.

  1. Bugs:
  • Indexing errors
  • Computes vert_vort_mag when using Leith+E
  • Fixes some logic in the smoothing code
  • Fixes application of upper bounds on Ah
  1. Features:
  • Smoothing subroutine (smooth_x9) has been re-written to preserve rotational symmetry.
  • Computations re-arranged to improve efficiency (fewer blocking halo updates)
  • A new parameter SMOOTH_AH added. Default is True. When False, no smoothing is applied to the backscatter and biharmonic coefficients, which significantly reduces the number of blocking halo updates.
  • Added some code to make Leith+E work with open boundaries.

iangrooms and others added 9 commits October 12, 2023 19:05
Biharmonic Leith uses Del omega at is-1 and js-1. This unavoidably requires
u at js-3 and v at is-3, which are unavailable. It also requires Del omega
at Ieq+1 and Jeq+1, which requires v at Ieq+3 and u at Jeq+3, which are
unavailable. This necessitates a halo update.

Fixes several bugs in Leith+E.
- Fixes indexing when computing smoothed vorticity and its gradient
- Crucially, computes `vert_vort_mag` when using Leith+E
- Fixes some logic in the smoothing code
- Other minor indexing fixes
Ah is required at h and q points. The original code computed Ah at
h points, then packed into Ah_h, then applied upper bounds to Ah.
If Ah_h is in the diag_table or if debug is true, then the value of
Ah with upper bounds get packed into Ah_h. Then, at q points the
code unpacks Ah_h. This update makes sure that the upper bound
gets applied to q points, not just h points.
The main thing that this commit does is to perform smoothing of u and v
outside of the loop over layers. This swaps nz 2D blocking halo updates
for a single blocking 3D halo update.
This commit adds a runtime flag, SMOOTH_AH. If True (default) then
`m_leithy` and `Ah` are both smoothed, which leads to many blocking
communications. If False then these fields are rougher, but there
is less communication.
This commit removes one halo update in Leith+E. To achieve this
requires re-indexing two assignments. The value of Ah and Kh are
computed at h points, then re-used at q points. Without the halo
update it is necessary to offset the assignment at h and q
points, e.g. Kh(I,J) = Kh_h(i+1,j+1,k), to avoid accessing
values that have not been computed.
Adds code so that Leith+E works with OBC.
This commit eliminates one more halo update in Leith+E.
  This commit revises the smoothing code used when USE_LEITHY = True to give
answers that respect rotational symmetry and it also corrects some horizontal
indexing bugs and problems with the staggering in some halo update and smooth_x9
calls and reduces some loop ranges to their minimal required values.  The
specific changes include:

  1. Corrected a horizontal indexing bug when interpolating Kh_h and Ah_h to
     corner (q) points when USE_LEITHY = True.  These had previously been
     inappropriately copied from the thickness point to the southwest of the
     corner point.  This required symmetric-memory-mode calculations of the
     thickness point viscosities whenever USE_LEITHY is true, but to avoid adding
     complicated logic, the symmetric-memory loop bounds are used for the
     calculation of Kh.

  2. Revised smooth_x9 to give rotationally symmetric answers and split it into
     the two routines smooth_x9_h and smooth_x9_uv to reduce the memory used by
     this routine and reduce the use of optional arguments.

  3. Eliminated 4 unneeded halo update calls, and added error handling for the
     case where Leith options are used with insufficiently wide halos.

  4. Added new integers to indicate the loop ranges over which the viscosities
     and related variables should be calculated, depending on which options are
     active, and then adjusted 91 do-loop extents horizontal_viscosity code to
     reflect the loop ranges over which arrays are actually used.

  5. Added a new 2-d variable for the squared viscosity for smoothing that can
     be used for halo updates and to avoid having a variable with confusingly
     inconsistent dimensions at various points in the code.

  6. Corrected the position arguments on 2 smooth_x9 calls and 4 pass_var calls
     that are used when USE_LEITHY=.true. and SMOOTH_AH=.true.  As previously
     written, these smooth_x9 and pass_var calls would work when in non-symmetric
     memory mode but would give incorrect answers when in symmetric memory mode.

  These revisions change answers when USE_LEITHY is true, but answers are
bitwise identical in all other cases.
+(*)Correct rotational symmetry with USE_LEITHY
@gustavo-marques
Copy link
Collaborator

Testing: pr_mom.derecho. No answer changes.

@gustavo-marques gustavo-marques self-requested a review February 10, 2024 23:03
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.

3 participants