Change shear mixing to LMD94 formulation and add ri smoothing by katsmith133 · Pull Request #352 · E3SM-Project/Omega

katsmith133 · 2026-02-23T16:59:35Z

This PR changes the shear vertical mixing formulation from Pacanowski and Philander (1981) to Large et al (1994), adds in vertical smoothing to the gradient Richardson number, and updates the tests and docs to reflect these changes.

Passes all CTests on PM-CPU and PM-GPU and the polaris omega_pr suite on PM.

Checklist

Testing

CTest unit tests

Machine: PM-CPU, PM-GPU
Compiler: gnu, gnugpu
Build type: Relase
Result: All tests passed

Polaris `omega_pr` suite

Baseline workdir: /pscratch/sd/k/katsmith/polaris_testing_vert_mixing//baseline_omega_pr
Baseline build: /pscratch/sd/k/katsmith/polaris-main/e3sm_submodules/Omega
PR build: /pscratch/sd/k/katsmith/polaris-pr352/e3sm_submodules/Omega
PR workdir: /pscratch/sd/k/katsmith/polaris_testing_vert_mixing/pr352_omega_pr
Machine: pm-cpu
Compiler: gnu
Build type: Release
Log: not found
Result: All tests passed

Tests added/modified/impacted

CTest: in VertMixTest added tests: testGradRichNum and testOneTwoOneFilter, and modified: testShearVertMix and testTotalVertMix.

Documentation

Documentation has be built locally here and charges are as expected: https://portal.nersc.gov/project/e3sm/katsmith/omega_docs/html/devGuide/VerticalMixingCoeff.html

vanroekel

Looks pretty good @katsmith133 just a few minor comments in the docs and one larger one in the main code.

katsmith133 · 2026-02-25T16:52:32Z

Thanks for the review, @vanroekel! I will have a look at it today. Seems like I misspelled LMD94 A LOT, haha.

katsmith133 · 2026-02-26T17:20:38Z

@cbegeman I mostly included you on this review because the polaris omega_pr tests just hang and I want to figure out if this is a polaris problem or a this PR problem, especially after Youngsung mentioned that he saw these tests hanging on PM as well. Let me know your thoughts.

cbegeman · 2026-02-26T17:52:00Z

@katsmith133 Thanks for the ping. I'd like to do polaris tests including the single column shear mixing one. Can you provide me with the appropriate mapping from Omega config options to MPAS-O config options (or indicate where there is no direct mapping? See https://github.com/E3SM-Project/polaris/pull/435/files#diff-d953af580e6050fd590b94c370ebc3401d6ed1b24e282d8d2d568536c35f8949R132 for an example.

Is there anything else I should know about when comparing to MPAS-O and with the current Omega develop?

alicebarthel · 2026-06-04T17:11:23Z

-dependent shear instability driven mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These are linearly additive and are describe a bit
-more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future.
+background mixing value, (2) a convective instability mixing value, and (3) a Richardson number-dependent shear-instability-driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 shear instability driven mixing parameterization. These are linearly additive and are describe a bit
+more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear instability driven mixing calculation.


Suggested change

more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear instability driven mixing calculation.

more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 (vertical) filter before using it in the shear instability driven mixing calculation.

alicebarthel · 2026-06-04T18:09:42Z

+                Kokkos::max(0.0_Real,
+                            0.5_Real * (BruntVaisalaFreqSq(ICell, K2) +
+                                        BruntVaisalaFreqSq(JCell, K2))) /
+                (ShearSquared + 1.0e-12_Real);


@mwarusz Is it better practice to do:
(ShearSquared + 1.0e-12_Real)
or Kokkos::max(ShearSquared, Eps) ?
I had that question for other parts of my code under development.

I am not sure this can be answered in general . (ShearSquared + 1.0e-12_Real) introduces large errors for values of ShearSquared on the order of 1e-12. Does this matter in this expression ? Kokkos::max(ShearSquared, Eps) doesn't have this problem, but it could be slower. Hard to be sure without benchmarking though.

This seems like a design choice much larger than this PR.

I would suggest not worrying about the difference between the two, but I'd wonder if we want to have hard coded values like this. Maybe instead we should have (ShearSquared + Eps) and define Eps = 1.0e-12_Real elsewhere?

@vanroekel do you think Eps should be defined in a central place for all of Omega to use (maybe GlobalConstants.h)? Not sure if something like 1.0e-12_Real is used elsewhere, but I could do a quick PR that replaces all of them with Eps

Yes I'd prefer we have a common small value for denominators that is unified if possible

@vanroekel Ok, thanks! I've created an issue for this and will address it with a separate PR that touches all instances of needing a small value (I don't think there are many yet).

alicebarthel · 2026-06-10T14:43:51Z

Just a summary of my review.
The 2 themes that came up go beyond the scope of this particular PR, and we should probably agree on how we proceed:

Do we require new PRs to use minLayerCell or are we ok with them assuming 0? We originally let code assume 0 since the non-zero case will be when we introduce sub-ice-shelf cavities. If this one assumes 0, it will need to be revised then. So do we ask for the adjustment now or later?
What is the fill values? This goes beyond vertMix and is being discussed/addressed in Fix fill values in all Omega fields #428.

Apart from these 2 points, my review is pretty minor. I don't feel very confident about having checked the edge cases for the loop limits but maybe @mwarusz has gone through carefully?

cbegeman · 2026-06-10T15:24:13Z

Regarding #352 (comment) and

1. Do we require new PRs to use minLayerCell or are we ok with them assuming 0? We originally let code assume 0 since the non-zero case will be when we introduce sub-ice-shelf cavities. If this one assumes 0, it will need to be revised then. So do we ask for the adjustment now or later?

I just want to contribute my experience with removing the assumption of KMin = 0 in MPAS-O. It can be pretty difficult to locate all the places this assumption is made after the fact so I really support using the minLayerCell variable now.

katsmith133 · 2026-06-10T15:37:35Z

Regarding #352 (comment) and
1. Do we require new PRs to use minLayerCell or are we ok with them assuming 0? We originally let code assume 0 since the non-zero case will be when we introduce sub-ice-shelf cavities. If this one assumes 0, it will need to be revised then. So do we ask for the adjustment now or later?
I just want to contribute my experience with removing the assumption of KMin = 0 in MPAS-O. It can be pretty difficult to locate all the places this assumption is made after the fact so I really support using the minLayerCell variable now.

@cbegeman I think this makes sense, thanks for your thoughts on this! Given that we have several places this change already needs to be made, do you think it would be better to have one big PR (like the PseudoThickness PR) that goes through and makes all of the current changes in one go? Or just rely on each new PR chipping away at the issue over time? I feel like the first would be better to get it all, but also takes someone's time away that is a rather precious resource right now. Perhaps a discussion for the new group meeting?

@alicebarthel and @mwarusz, thanks for the additional review comments! I will get started on making the changes and responding to them shortly.

sbrus89

CTests and omega_pr suite passed (with baselines) with a local merge on Frontier (craygnu).

Test project /ccs/home/brus/run/omega_ctest/build_omega/build_frontier_craygnu
      Start  1: DATA_TYPES_TEST
 1/40 Test  #1: DATA_TYPES_TEST ............................   Passed    2.35 sec
      Start  2: MACHINE_ENV_TEST
 2/40 Test  #2: MACHINE_ENV_TEST ...........................   Passed    0.29 sec
      Start  3: BROADCAST_TEST
 3/40 Test  #3: BROADCAST_TEST .............................   Passed    0.28 sec
      Start  4: LOGGING_TEST
 4/40 Test  #4: LOGGING_TEST ...............................   Passed    0.28 sec
      Start  5: ERROR_TEST
 5/40 Test  #5: ERROR_TEST .................................   Passed    0.30 sec
      Start  6: DECOMP_NTASK1_TEST
 6/40 Test  #6: DECOMP_NTASK1_TEST .........................   Passed    0.61 sec
      Start  7: DECOMP_NTASK8_TEST
 7/40 Test  #7: DECOMP_NTASK8_TEST .........................   Passed    0.52 sec
      Start  8: HALO_TEST
 8/40 Test  #8: HALO_TEST ..................................   Passed    0.62 sec
      Start  9: DIMENSION_TEST
 9/40 Test  #9: DIMENSION_TEST .............................   Passed    0.53 sec
      Start 10: HORZMESH_TEST
10/40 Test #10: HORZMESH_TEST ..............................   Passed    1.19 sec
      Start 11: HORZOPERATORS_PLANE_TEST
11/40 Test #11: HORZOPERATORS_PLANE_TEST ...................   Passed    1.12 sec
      Start 12: HORZOPERATORS_SPHERE_TEST
12/40 Test #12: HORZOPERATORS_SPHERE_TEST ..................   Passed    1.23 sec
      Start 13: AUXVARS_PLANE_TEST
13/40 Test #13: AUXVARS_PLANE_TEST .........................   Passed    1.11 sec
      Start 14: AUXVARS_SPHERE_TEST
14/40 Test #14: AUXVARS_SPHERE_TEST ........................   Passed    1.14 sec
      Start 15: AUXSTATE_TEST
15/40 Test #15: AUXSTATE_TEST ..............................   Passed    1.29 sec
      Start 16: IO_TEST
16/40 Test #16: IO_TEST ....................................   Passed    1.95 sec
      Start 17: CONFIG_TEST
17/40 Test #17: CONFIG_TEST ................................   Passed    0.30 sec
      Start 18: FIELD_TEST
18/40 Test #18: FIELD_TEST .................................   Passed    0.53 sec
      Start 19: IOSTREAM_TEST
19/40 Test #19: IOSTREAM_TEST ..............................   Passed   14.62 sec
      Start 20: TEND_PLANE_TEST
20/40 Test #20: TEND_PLANE_TEST ............................   Passed    2.19 sec
      Start 21: TEND_PLANE_SINGLE_PRECISION_TEST
21/40 Test #21: TEND_PLANE_SINGLE_PRECISION_TEST ...........   Passed    2.24 sec
      Start 22: TEND_SPHERE_TEST
22/40 Test #22: TEND_SPHERE_TEST ...........................   Passed    2.22 sec
      Start 23: TENDENCIES_TEST
23/40 Test #23: TENDENCIES_TEST ............................   Passed    2.23 sec
      Start 24: PGRAD_TEST
24/40 Test #24: PGRAD_TEST .................................   Passed    2.23 sec
      Start 25: STATE_TEST
25/40 Test #25: STATE_TEST .................................   Passed    1.41 sec
      Start 26: TIMEMGR_TEST
26/40 Test #26: TIMEMGR_TEST ...............................   Passed    0.25 sec
      Start 27: TIMEINTERVAL_PARSE_TEST
27/40 Test #27: TIMEINTERVAL_PARSE_TEST ....................   Passed    0.15 sec
      Start 28: TIMEINTERVAL_PARSE_EXTENDED_FORMATS_TEST
28/40 Test #28: TIMEINTERVAL_PARSE_EXTENDED_FORMATS_TEST ...   Passed    0.13 sec
      Start 29: REDUCTIONS_TEST
29/40 Test #29: REDUCTIONS_TEST ............................   Passed    0.32 sec
      Start 30: TIMESTEPPER_TEST
30/40 Test #30: TIMESTEPPER_TEST ...........................   Passed    1.29 sec
      Start 31: KOKKOS_FLATPAR_TEST
31/40 Test #31: KOKKOS_FLATPAR_TEST ........................   Passed    0.21 sec
      Start 32: KOKKOS_HIPAR_TEST
32/40 Test #32: KOKKOS_HIPAR_TEST ..........................   Passed    0.16 sec
      Start 33: DRIVER_TEST
33/40 Test #33: DRIVER_TEST ................................   Passed    2.00 sec
      Start 34: TRACERS_TEST
34/40 Test #34: TRACERS_TEST ...............................   Passed    1.12 sec
      Start 35: TRIDIAGSOLVERS_TEST
35/40 Test #35: TRIDIAGSOLVERS_TEST ........................   Passed    0.27 sec
      Start 36: GSWC_CALL_TEST
36/40 Test #36: GSWC_CALL_TEST .............................   Passed    0.03 sec
      Start 37: EOS_TEST
37/40 Test #37: EOS_TEST ...................................   Passed    1.12 sec
      Start 38: VERTCOORD_TEST
38/40 Test #38: VERTCOORD_TEST .............................   Passed    1.20 sec
      Start 39: VERTMIX_TEST
39/40 Test #39: VERTMIX_TEST ...............................   Passed    1.10 sec
      Start 40: VERTADV_TEST
40/40 Test #40: VERTADV_TEST ...............................   Passed    1.07 sec

100% tests passed, 0 tests failed out of 40

Label Time Summary:
Omega-0    =  53.18 sec*proc (39 tests)
SERIAL     =  53.18 sec*proc (39 tests)

Total Test time (real) =  69.67 sec

Polaris `omega_pr` suite

Baseline workdir: /ccs/home/brus/run/polaris_omega_pr_baseline
Baseline build: /ccs/home/brus/run/polaris_omega_pr_baseline/build
PR build: /ccs/home/brus/run/polaris_vert_mix_omega_pr/build
PR workdir: /ccs/home/brus/run/polaris_vert_mix_omega_pr
Machine: frontier
Partition: batch
Compiler: craygnu
Build type: Release
Log: not found
Result: All tests passed

Thanks for all your work on this @katsmith133! This seems ready to merge unless you have additional changes to push or @mwarusz or @alicebarthel have further comments.

mwarusz · 2026-06-17T22:30:55Z

+               ComputeOneTwoOneFilter); /// Local view for 1-2-1 filter
   OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
   OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);



Suggested change

OMEGA_SCOPE(LocBackDiff, BackDiff);

OMEGA_SCOPE(LocBackVisc, BackVisc);

mwarusz · 2026-06-17T22:31:20Z

+                    LocVertDiff(ICell, K) = BackDiff;
+                    LocVertVisc(ICell, K) = BackVisc;


Suggested change

LocVertDiff(ICell, K) = BackDiff;

LocVertVisc(ICell, K) = BackVisc;

LocVertDiff(ICell, K) = LocBackDiff;

LocVertVisc(ICell, K) = LocBackVisc;

alicebarthel · 2026-06-17T23:16:56Z

I did a quick read-through and that looks good to me. Thanks @katsmith133 for your work on this! and to every one else for your reviews, comments, and general engagement.

katsmith133 requested review from alicebarthel, cbegeman and vanroekel February 23, 2026 16:59