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

Some touch-ups to porous barrier #186

Closed
wants to merge 13 commits into from

Conversation

herrwang0
Copy link

@herrwang0 herrwang0 commented Aug 7, 2022

This PR introduces some small changes to the porous barrier feature.
Edit: There are a number of small commits, maybe a squash commit should be done instead?

New features:

  • Added a new subroutine in MOM_fixed_initialization to read global sub-grid scale topography parameters at velocity points.
  • Added a new parameter to choose different interpolation method for velocity point interface height.
  • Diagnostics of the weight functions, CPU clocks and chksum debugs are added.
  • A control structure is added in module MOM_porous_barriers to manage relevant input parameters, diagnostics, and etc.
  • Added a global switch for porous barrier functionality to prevent unnecessary calculations and stacks when this feature is not used. The default at the moment is .True. (switched on) though, for backward compatibility.

Bug fixes:

  • Previously there is a bug in loop range in which the velocity points on the north/east edges of the domains are not calculated and there was not halo update to cover these points. This is fixed by change the loop range from data domain to computation domain + halo update
  • There was a bug with potentially incorrect eta_prev and A_layer_prev values, which was caused by a misplace of their updates inside an if-block. This is fixed after restructuring the corresponding code blocks
  • There was an implicit assumption that porous barrier is only applied to ocean depth below 0.0, which precluded wetting/drying. This is fixed by adding an additional parameter PORBAR_MASKING_DEPTH to control this critical depth.

Refactors:

  • Grouped porous barrier variables in MOM_variables is change from pointer to allocatable.
  • The calculations of the interface and layer weight functions are now split. Currently, there are two update of porous barriers in the step_MOM, which in fact target the layers and interfaces separately. In the future, it might be helpful to use fortran interface to unify the subroutine names.
  • The loops are re-arranged to increase readability.
  • An updating do_next mask is added to skip weight function calculation for all points above the shallowest point. This can make multi-layer runs (most layers in the open ocean are shallower than the topography) a bit more efficient.
  • An upper bound of 1.0 is applied to the layer-averaged weight function.
  • The above two points do give answers that are not bit-for-bit with the old runs. For consistency purposes, the newly-introduced ANSWER_DATE flag is used to recover old results. (Question: the current value of PORBAR_ANSWER_DATE is 20220806, should it be instead some future date when this PR is merged? )

@codecov
Copy link

codecov bot commented Aug 7, 2022

Codecov Report

Merging #186 (b296aa5) into dev/gfdl (e3c71b0) will increase coverage by 0.04%.
The diff coverage is 32.87%.

❗ Current head b296aa5 differs from pull request most recent head 97b00c7. Consider uploading reports for the commit 97b00c7 to get more accurate results

@@             Coverage Diff              @@
##           dev/gfdl     #186      +/-   ##
============================================
+ Coverage     37.22%   37.26%   +0.04%     
============================================
  Files           261      261              
  Lines         72340    72409      +69     
  Branches      13537    13532       -5     
============================================
+ Hits          26929    26984      +55     
- Misses        40417    40425       +8     
- Partials       4994     5000       +6     
Impacted Files Coverage Δ
src/core/MOM_CoriolisAdv.F90 50.74% <ø> (ø)
src/core/MOM_continuity.F90 60.71% <ø> (ø)
src/core/MOM_continuity_PPM.F90 72.54% <ø> (ø)
src/core/MOM_dynamics_split_RK2.F90 65.56% <ø> (ø)
src/core/MOM_dynamics_unsplit.F90 79.43% <ø> (ø)
src/core/MOM_dynamics_unsplit_RK2.F90 81.06% <ø> (ø)
src/core/MOM_grid.F90 70.89% <ø> (ø)
src/core/MOM_variables.F90 46.96% <ø> (ø)
src/framework/MOM_dyn_horgrid.F90 55.66% <ø> (ø)
src/initialization/MOM_shared_initialization.F90 29.29% <0.00%> (-1.25%) ⬇️
... and 8 more

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@herrwang0 herrwang0 closed this Aug 7, 2022
@herrwang0 herrwang0 reopened this Aug 8, 2022
Copy link
Member

@Hallberg-NOAA Hallberg-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that this PR adds some useful new capabilities and efficiencies for MOM6.
I have visually reviewed all of these changes, and they all look good to me. Although there are 12 substantive commits in this PR, each is a useful contribution in its own right, so they might we worth keeping, especially if this whole PR can synchronized with dev/gfdl via a rebase. Otherwise, some of the individual commits date back several months, so keeping the individual commits might not make sense for debugging any issues, in which case this PR should be dealt with via a squash merge.

I am approving this commit based on this visual review, conditional upon the updated TC testing and pipeline testing passing.

I would, however, be interested in any comments @sditkovsky might have on this PR.

@sditkovsky
Copy link

@Hallberg-NOAA Thanks for the notification. I have rerun some of my test cases with these changes and the results look good to me. And I think it's great to have the new options for read-in method and interpolation schemes.

@herrwang0 Thanks for the work. Glad the porous barriers functionality is coming in handy!

@marshallward
Copy link
Member

@herrwang0 This is ready to be merged, but there are conflicts in some of the commits and it cannot be rebased.

If you can be troubled to do a rebase (git rebase -i dev/gfdl), fix the conflicts, and then force-push an update, that would be appreciated. But if it's too much trouble, then we can just also do a regular (FF) merge.

Completely up to you, no pressure.

* Porous barrier variables, which are grouped in porous_barrier_ptrs,
are changed from pointers to allocatables.
* Copies (aliases) of these variables in MOM_CS are removed.
* The name porous_barrier_ptrs is yet to be changed to minimize the
number of files needed to be edited.
* The interface porous_widths is deleted and subroutine por_widths is
renamed as porous_widths.
* porous_barrier_CS is added to control input and diagnostics of the
 module
* Diagnostics for both interface and layer averages weights are added in
subroutine porous_widths.
* An _init subroutine is added to facilitate reading parameters and
registering diagnostics
* checksum debugs are added within subroutine porous_widths.
This commit primarily fixes indexing bugs in subroutine porous_widths.

* Loop range is subroutine porous_widths is changed from data domain to
computation domain.
* find_eta call is now with a proper halo to cover all velocity points.
* Halo update for porous barrier fields is added in step_MOM_*.

Other changes:
* mask_depth, a component of porous_barrier_CS is now used to as
the criterion for masking. This removes the limit that the cell edge
depth has to be below sea surface.
* Output variable eta_cor was unused and is now changed to a local.
* Some unused indexing variables are removed.
parameters

* subroutine set_subgrid_topo_at_vel_from_file is added to read max, min
and avg depth at the edges from file.
* The subroutine is only called when a new runtime parameter
SUBGRID_TOPO_AT_VEL is True. Default is false.
* id_clock is added (as a private variable) to time porous barrier
calculations.
* Some small format fixes
* A runtime parameter PORBAR_ETA_INTERP is added to control different
methods for calculating interface height at the velocity points from
adjacent tracer cells.
* Two small thickness variables are added and scaled the unit of eta.
This is to assit the harmonic mean calculation, but eventually the
if-test eta_s - eta_prev>0 should be relaced by Angstrom.
* Runtime parameter POROUS_BARRIER_MASKING_DEPTH is renamed to make it
a liitlle bit shorter.
* log_version call is added to separate out porous_barriers module in
parameter file.
The main loops in porous_widths are simpified and cleaned up.
* calc_por_layer iss slightly modified to reduced the if-blocks in
the main loops.
* k-loops are moved out.
* To assist this new structure, two 2-D arrays are added to the stack.
* A new function eta_at_uv is added to treat different interpolations.
This is currently done at every point within the loop, and it is
probably adding too much an overhead. It is better pre-calculate
 eta_[uv] all together, which would increase the stack size.
The averaged weight over a layer (por_face_area[UV]) and the weight
at the interface (por_layer_width[UV]) are now calculated
separately, as the only two updates in step_mom are in fact
using only one of them. This reduces
some unnecessary calculations.

* Two new subroutines replaced the original `porous_widths`. They could
be combined with interface if we remove porous_barrier_ptrs type, and
use variables (of different nz) as output.
* There is a place holder switch do_next in the two calc_por subs. This
is used to further reduce calculations for all layers above D_max.
Will be tested and implemented in the next commit.
* A new subroutine calc_eta_at_uv is added to calculate eta at uv.
This simplifies the code, but also increases stack and some performance
overhead.
* A new runtime parameter USE_POROUS_BARRIER is added to control this
feature. The default is True at the moment to assist potential test. It
should be changed to false in the future.
* A boolean array is added to avoid unnecessary porous barrier
 calculations above the shallowest points.
* The layer-averaged weights are bounded by 1.0 explicitly, to avoid
the cases it goes beyond due to roundoff errors.
* For very thin layers (Angstrom), layer-averaged weights are set to
zeros for simplicity.
  * Perhaps it is more reasonable to use the inferace weight for these
  cases.
  * It might be useful to make the thin layer definition a runtime
  parameter.
* A bug related the size of eta_[uv] is fixed.

This commit is potentially answer-changing when the porouss barrier
module is used.
* The recently introduced ANSWER_DATE functionality is used to maintain a
version that should be bit-for-bit reproducible for previous porous barrier
related runs.
* Rename porous_barrier_ptrs to porous_barrier_type
* Some small documentation edits
* Fixed a bug that eta_[uv] was not properly initialized, which caused issues
with global chksums with DEBUG=True.
* Documents added to comply with Doxygen tests
@herrwang0 herrwang0 force-pushed the porous_barrier_touchup branch from 89f239e to 9848f17 Compare August 8, 2022 20:50
@marshallward
Copy link
Member

This appears to be failing in our static-memory builds in our CI. This is the backtrace:

#1  0x955cf8 in __mom_MOD_initialize_mom
	at ../../../../../.././/src/MOM6/src/core/MOM.F90:2490
#2  0xc34718 in mom6
	at ../../../../../.././/src/MOM6/config_src/drivers/solo_driver/MOM_driver.F90:284
#3  0xc368fa in main
	at ../../../../../.././/src/MOM6/config_src/drivers/solo_driver/MOM_driver.F90:27

and the offending line is

ALLOC_(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0

with the following error:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

I haven't yet probed further, but it is possible that one of these indices is not correctly being assigned its static-build index value.

@marshallward
Copy link
Member

marshallward commented Aug 18, 2022

@herrwang0 I believe the issue is that the contents of porous_barrier_type were declared as dynamic arrays:

type, public :: porous_barrier_type
! Each of the following fields has nz layers.
real, allocatable :: por_face_areaU(:,:,:) !< fractional open area of U-faces [nondim]
real, allocatable :: por_face_areaV(:,:,:) !< fractional open area of V-faces [nondim]
! Each of the following fields is found at nz+1 interfaces.
real, allocatable :: por_layer_widthU(:,:,:) !< fractional open width of U-faces [nondim]
real, allocatable :: por_layer_widthV(:,:,:) !< fractional open width of V-faces [nondim]
end type porous_barrier_type

but that you are trying to allocate them with the static memory macro ALLOC_ which does nothing in static mode (since it assumes that the array is already declared at fixed size.)

MOM6/src/core/MOM.F90

Lines 2490 to 2493 in a14c292

ALLOC_(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0
ALLOC_(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%pbv%por_face_areaV(:,:,:) = 1.0
ALLOC_(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%pbv%por_layer_widthU(:,:,:) = 1.0
ALLOC_(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%pbv%por_layer_widthV(:,:,:) = 1.0

You will need to either use dynamic memory syntax for these fields, or re-declare the fields with the static memory macros (real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: ...).

(My preference is always for dynamic memory, but either is suitable.)

@herrwang0
Copy link
Author

Thanks @marshallward for catching this!

The mix-up really exposes that I had no idea which was preferable!

My understanding is that ALLOC_, ALLOCABLE_ ,NIMEM_, and etc defined in the memory macro cover both static and dynamic memory, is that right? Does that mean we should use this set?
But I suppose it would prohibit the use of features like allocate(..., source= )

@adcroft
Copy link
Member

adcroft commented Aug 19, 2022

Compare your declarations to that of CS%ssh_rint. You'll see that you used the _PTR_ macro, meant for pointers, Easy oversight - I do it all the time!

@herrwang0 herrwang0 force-pushed the porous_barrier_touchup branch from a14c292 to bbc59f8 Compare August 22, 2022 18:25
@herrwang0
Copy link
Author

I decided to simply change it to the dynamic memory allocate. It might be more beneficial to switch these variables to subroutine input/output, but that could be a discussion in the future.

@herrwang0
Copy link
Author

Compare your declarations to that of CS%ssh_rint. You'll see that you used the _PTR_ macro, meant for pointers, Easy oversight - I do it all the time!

Hi @adcroft , you were probably looking at the previous version of it

  real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: por_face_areaU !< fractional open area of U-faces [nondim]
  real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: por_face_areaV !< fractional open area of V-faces [nondim]
  real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: por_layer_widthU !< fractional open width
                                                                                   !! of U-faces [nondim]
  real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: por_layer_widthV !< fractional open width
                                                                                   !! of V-faces [nondim]

In this PR, this part is actually now gone and the grouped type (defined in MOM_variables) is used instead.

Sorry it took this long to reply, I just realized what you were referring to just now (when I was solving the merge conflicts).

@marshallward
Copy link
Member

marshallward commented Aug 24, 2022

@herrwang0 I have a rebased version of this PR, could you review it here?

https://github.com/marshallward/MOM6/tree/pr_186_rebase

For comparison:
dev/gfdl...marshallward:MOM6:pr_186_rebase

(We also might want to combine some of these commits, such as the bugfixes)

@herrwang0
Copy link
Author

@marshallward Thanks for doing the rebasing (I should have done that instead of merge before the latest bug fix)!
It looks good to me. Yes, I agree that these bugfixes should be combined.

@marshallward
Copy link
Member

Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/16574 ✔️ 🟡

The rebased version has passed, so I will merge that one.

@marshallward
Copy link
Member

These commits have been merged in a rebased branch, so closing this PR.

Note that the commits will have different hashes.

@herrwang0 herrwang0 deleted the porous_barrier_touchup branch May 20, 2024 19:30
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.

5 participants