Implement P-star Vertical Coordinate by brian-oneill · Pull Request #401 · E3SM-Project/Omega

brian-oneill · 2026-05-05T21:34:55Z

Changes include:

Read NVertLayers from InitVertCoord stream
Rename RefLayerThickness to RefPseudoThickness to avoid conflation with MPAS field refLayerThickness
Add Field metadata for RefPseudoThickness to InitVertCoord stream
Needed updates to IOStream class
Add computeTargetThickness method to computeMomVertAux
Add ALE contribution to computeVerticalVelocity

Checklist

Fixes #391

brian-oneill · 2026-05-05T21:48:31Z

Testing

Machines: pm-cpu, pm-gpu, frontier
Compilers: gnu, gnugpu, craycray, craycray-mphipcc
Build type: Release
Result: All tests passed

Polaris `omega_pr` suite

Baseline workdir: /lustre/orion/proj-shared/cli115/boneill/polaris/testing/baseline_omega_pr
Baseline build: /lustre/orion/proj-shared/cli115/boneill/polaris/polaris/e3sm_submodules/Omega/_build
PR build: /lustre/orion/proj-shared/cli115/boneill/myFork/prCPU
PR workdir: /lustre/orion/proj-shared/cli115/boneill/polaris/testing/my_branch_omega_pr
Machine: frontier
Partition: batch
Compiler: craygnu
Build type: Release
Log: /lustre/orion/proj-shared/cli115/boneill/polaris/testing/my_branch_omega_pr/polaris_omega_pr.o4531781
Result: All tests passed

xylar · 2026-05-06T05:57:11Z

@brian-oneill, thanks for making these changes! I'll give this a look.

But, in response to your question, yes, we want to read in VertCoordMovementWeights as part of the vertical initial condition stream. See the vertical coordinate design doc:
https://docs.e3sm.org/Omega/Omega/design/VertCoord.html#creation
The design is missing MinLayerCell and MaxLayerCell but we have those already in the stream.

xylar

@brian-oneill, this looks excellent to me. I just have two comments based on reviewing the code itself, both quite minor. I will go on to revising my Polaris branch with your (excellent) variable name change and make sure I see passing tests, the same as you.

xylar · 2026-05-06T10:22:43Z

Excellent! I'm also seeing that the omega_pr suite passes with these changes and E3SM-Project/polaris#573:

Polaris `omega_pr` suite

Baseline workdir: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260504/omega-pr-baseline
Baseline build: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260504/omega-pr-baseline/build
PR build: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260506/omega-pr-p-star/build
PR workdir: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260506/omega-pr-p-star
Machine: frontier
Partition: batch
Compiler: craygnu
Build type: Release
Log: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260506/omega-pr-p-star/polaris_omega_pr.o4533646
Result: All tests passed

Once my two small concerns above get addressed, I'll be happy to approve and merge.

xylar · 2026-05-06T10:23:22Z

Oh, and after the docs get updated. Presumably the changes there won't be large.

xylar · 2026-05-06T11:47:02Z

@brian-oneill, when you update the docs, could you take out or replace the references to LayerThicknessPStar in the vertical coordinate docs? That variable doesn't seem to exist.

xylar · 2026-05-06T14:34:23Z

@hyungyukang, can you check if this works with the overflow test case as part of a review here?

xylar · 2026-05-06T14:34:58Z

@sbrus89, if it's feasible, we should test this with the baroclinic channel.

cbegeman · 2026-05-06T18:35:38Z

                    const I4 K = KStart + KVec;
                    LocLayerThickTarget(ICell, K) =
-                        LocRefLayerThick(ICell, K) *
+                        LocRefPseudoThick(ICell, K) *


Does this imply that LocLayerThickTarget is actually in pseudo-thickness, not geometric thickness? If so, should we rename it LocPseudoThickTarget and should AleTerm above be using (LocThickTarget - PseudoThickness) or convert LocLayerThickTarget to geometric thickness?

This is not the purpose of this PR. I will rename this in #327

If we rename things other than *RefLayer* to *RefPseudo*, there will be essentially no end in what needs fixing.

I took the initiative to rename RefLayer to RefPseudo specifically because the reference thickness is added to the initialization stream in this PR and I just felt uneasy keeping the name as is, since refLayerThickness exists as a distinct 1D field in the MPAS meshes we've been using for tests.

As for the ALE term, all these quantities are pseudo-thicknesses, and the vertical velocities we use in Omega are pseudo-velocities $\tilde{w} \equiv d\tilde{z} / dt$, so no conversion is necessary.

Thanks for clarifying. Good to know that it's all pseudothickness here.

Ah, sorry for missing the nuance here. Thanks for clarifying it @brian-oneill.

brian-oneill · 2026-05-07T19:23:44Z

But, in response to your question, yes, we want to read in VertCoordMovementWeights as part of the vertical initial condition stream.

@xylar (and anyone else): it looks like MPAS reads the movement weights from the mesh, but also the config_vert_coord_movement option can override the values read from the mesh. Do we want to preserve this behavior in Omega, allowing the weights to be switched via the MovementWeightType config option, or should I remove that option and we rely solely on the weights read from the stream?

xylar · 2026-05-07T19:31:47Z

@mark-petersen, you might be best equipped to answer @brian-oneill's question. My feeling would be that we typically want to change VertCoordMovementWeights in idealized tests and that it would be easier and more consistent with how we do testing if this is something Polaris handles, not something Omega overrides. Do you foresee any circumstances where an Omega override would be important?

cbegeman · 2026-05-07T21:11:09Z

@brian-oneill For what it's worth, I'm comfortable with you removing this option and only reading it in from the mesh.

Only call IO::destroyDecomp for distributed fields in readFieldData

xylar · 2026-05-08T11:44:02Z

@sbrus89 and @hyungyukang, any update on how your testing is going on this? We need to keep this moving.

cbegeman · 2026-05-08T14:29:08Z

@sbrus89 and @hyungyukang, any update on how your testing is going on this? We need to keep this moving.

@hyungyukang Let me know if you'd like me to test the overflow case.

hyungyukang · 2026-05-08T14:32:18Z

@xylar, @cbegeman , I'm working on it now and will post the results soon!

hyungyukang · 2026-05-08T15:18:37Z

Here's the overflow test result:

It looks the same as before. Is that expected, or should I update something in omega.yml?

sbrus89 · 2026-05-08T15:45:39Z

I'm just finalizing the baroclinic_channel branch in polaris. I'll report back soon.

xylar · 2026-05-08T16:59:05Z

It looks the same as before. Is that expected, or should I update something in omega.yml?

@hyungyukang, my understanding is that something should have changed. We are now calling:
https://github.com/brian-oneill/E3SM/blob/8e98997e4317f0d93787dc2528bd5c88b7666dbc/components/omega/src/ocn/AuxiliaryState.cpp#L102-L103
which we weren't before and we are now computing vertical advection accounting for coordinate movement:
https://github.com/brian-oneill/E3SM/blob/8e98997e4317f0d93787dc2528bd5c88b7666dbc/components/omega/src/ocn/VertAdv.cpp#L442-L458
which, again, we weren't before. Together, these should mean that the vertical coordinate and vertical advection behave differently from before.

I don't see anything in your overflow configuration (https://github.com/E3SM-Project/polaris/pull/572/changes#diff-3b61ea3e604ca93614415174a1dc7c8844506c2d3052e119f26a56ac301e37aa) to explain why you wouldn't be seeing any differences.

@brian-oneill, any guess why this isn't changing answers?

xylar · 2026-05-08T17:00:32Z

@hyungyukang, the differences might be small. Have you verified that you're not seeing changes, or just that things look similar by eye?

hyungyukang · 2026-05-08T17:34:37Z

@hyungyukang, the differences might be small. Have you verified that you're not seeing changes, or just that things look similar by eye?

@xylar , thanks a lot for the explanations. It was just an eyeball comaprison. I will check the field differences of the overflow test between p* and z* and post here later.

xylar · 2026-05-08T17:37:43Z

between p* and z* and post here later

How are you getting z*? I don't think that's what Omega produces before this change, but I could be mistaken.

Or do you mean z* from MPAS-Ocean?

brian-oneill · 2026-05-08T20:26:35Z

Some meshes like the QU240 mesh (including both the old one being used for ctests and the new ocean.QU.240km.omega_vars.260506.nc mesh with the newer omega variable names) have all movement weights equal 0. This leads to division by 0 in computeTargetThickness.

This doesn't seem intentional here, but I'm wondering what the desired behavior should be when this is encountered. Should a mesh with all zero movement weights be treated as a flag for an impermeable interfaces configuration and disable vertical transport? Or should initialization default to p* when invalid weights like this are read? Or just abort with an error?

cbegeman · 2026-05-08T21:46:37Z

Some meshes like the QU240 mesh (including both the old one being used for ctests and the new ocean.QU.240km.omega_vars.260506.nc mesh with the newer omega variable names) have all movement weights equal 0. This leads to division by 0 in computeTargetThickness.

This doesn't seem intentional here, but I'm wondering what the desired behavior should be when this is encountered. Should a mesh with all zero movement weights be treated as a flag for an impermeable interfaces configuration and disable vertical transport? Or should initialization default to p* when invalid weights like this are read? Or just abort with an error?

I can't speak to what we want for the QU mesh, but I think we do want to have a way to replicate the impermeable_interfaces behavior. There are many test cases that take advantage of this and I think it could be helpful for debugging. Without taking a close look at the code, it does seem like a zero-array movement weights would be a way to achieve this.

xylar · 2026-05-08T22:28:35Z

That seems like a glitch in the ocean.QU.240km.omega_vars.260506.nc and we want it to be all ones. I can fix that.

I think it should be an error if we have selected the p* coordinate and have all zero VertCoordMovementWeights. I think there should be a different way to specify impermeable layers that doesn't rely on p*, given that it seems like all zeros in VertCoordMovementWeights doesn't naturally lead to impermeable layers.

hyungyukang

@xylar , I meant the difference between this PR and develop. I was confusing it with MPAS-O’s z*.

Here is the difference from the 6-hour overflow test:

There are some differences, but they are small, as expected.

xylar

I'm approving based on a second visual inspection of the code and @hyungyukang's overflow testing.

xylar · 2026-05-11T12:04:05Z

+                           StreamName);
+            }
+            if (NumNonZero == 0) {
+               // TODO: ABORT_ERROR when all weights equal 0


Let's make an issue covering this once this PR gets merged. Let's make a second issue about needing a config option to specify impermeable layers.

xylar · 2026-05-11T12:06:22Z

                    const I4 K = KStart + KVec;
                    LocLayerThickTarget(ICell, K) =
-                        LocRefLayerThick(ICell, K) *
+                        LocRefPseudoThick(ICell, K) *


Ah, sorry for missing the nuance here. Thanks for clarifying it @brian-oneill.

xylar · 2026-05-11T12:07:36Z

+   // Projection factors for each RK stage. These are the RKC coefficients
+   // shifted by one RK stage and represent the fraction of the timestep over
+   // which fields are projected from the old state to the state at the end of a
+   // RK stage.


Thanks, this is very helpful!

xylar · 2026-05-11T12:10:05Z

+    TimeInterval ProjDt ///< [in] Time interval for projection over the current
+                        ///< time stepper stage


Thanks, that's helpful!

mark-petersen · 2026-05-11T20:49:00Z

My feeling would be that we typically want to change VertCoordMovementWeights in idealized tests and that it would be easier and more consistent with how we do testing if this is something Polaris handles, not something Omega overrides. Do you foresee any circumstances where an Omega override would be important?

I agree that VertCoordMovementWeights can be read on init, and Polaris can handle the generation. If the variable does not exist in the input file, Omega could issue a warning and default to all ones.

xylar · 2026-05-11T21:01:28Z

If the variable does not exist in the input file, Omega could issue a warning and default to all ones.

@mark-petersen, I think we agreed above that we would default to all ones temporarily but once files (e.g. those for CTests) get cleaned up, not having VertCoordMovementWeights would become an error, not just a warning.

cbegeman

Approving on the basis of visual inspection and the testing of others. @brian-oneill Thanks for your work on this!

brian-oneill added 13 commits May 4, 2026 13:03

Add getFilename function to IOStream

6a875e8

Read NVertLayers from InitVertCoord stream

dfa11af

Create calendar if needed in IOStream::init(void)

188d104

Add IOStream::init to tests where needed

a996d36

Rename RefLayerThickness to RefPseudoThickness

d09370c

Add RefPseudoThickness field

479f153

Exchange halo for RefPseudoThickness

87f9049

Add ALE contribution to vertical transport velocity

307c363

Add computeTargetThickness to computeMomVertAux

7c7b289

Propagate ProjDt through call tree to computeVerticalVelocity

4e4f39f

Add ALE projection time intervals to time steppers

0281bbd

Update ctests

ee31325

Correct TotalVerticalVelocity calculation

8e98997

brian-oneill mentioned this pull request May 6, 2026

VertCoord getting NVertLayers from wrong stream and not reading RefLayerThickness #391

Open

xylar reviewed May 6, 2026

View reviewed changes

Comment thread components/omega/src/infra/IOStream.cpp

Comment thread components/omega/src/ocn/Tendencies.cpp Outdated

xylar self-assigned this May 6, 2026

xylar requested a review from hyungyukang May 6, 2026 14:33

xylar requested a review from sbrus89 May 6, 2026 14:34

xylar reviewed May 6, 2026

View reviewed changes

Comment thread components/omega/src/timeStepping/RungeKutta4Stepper.cpp

cbegeman reviewed May 6, 2026

View reviewed changes

IOStream bug fix

81ba447

Only call IO::destroyDecomp for distributed fields in readFieldData

brian-oneill added 2 commits May 8, 2026 14:14

Add movement weights field to InitVertCoord stream

5830446

Remove MovementWeightType config option

096cdcf

sbrus89 mentioned this pull request May 8, 2026

Add omega support for baroclinic channel case E3SM-Project/polaris#576

Open

6 tasks

hyungyukang reviewed May 10, 2026

View reviewed changes

brian-oneill added 5 commits May 10, 2026 20:52

Check if all VertCoordMovementWeights are 0

c08d64d

Rename ProjCoeff

db60748

Add comment explaining RKProj coefficients

508c7e6

Update docs

9bdc244

Remove MovementWeightType option from default config

2b4b676

xylar approved these changes May 11, 2026

View reviewed changes

cbegeman approved these changes May 13, 2026

View reviewed changes

xylar mentioned this pull request May 13, 2026

Include Omega fields needed for p* coordinate in initial conditions E3SM-Project/polaris#573

Draft

1 task

		TimeInterval ProjDt ///< [in] Time interval for projection over the current
		///< time stepper stage

Conversation

brian-oneill commented May 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

brian-oneill commented May 5, 2026

Testing

Polaris omega_pr suite

Uh oh!

xylar commented May 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

xylar commented May 6, 2026

Polaris omega_pr suite

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar commented May 6, 2026

Uh oh!

Uh oh!

cbegeman May 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brian-oneill commented May 7, 2026

Uh oh!

xylar commented May 7, 2026

Uh oh!

cbegeman commented May 7, 2026

Uh oh!

xylar commented May 8, 2026

Uh oh!

cbegeman commented May 8, 2026

Uh oh!

hyungyukang commented May 8, 2026

Uh oh!

hyungyukang commented May 8, 2026

Uh oh!

sbrus89 commented May 8, 2026

Uh oh!

xylar commented May 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

xylar commented May 8, 2026

Uh oh!

hyungyukang commented May 8, 2026

Uh oh!

xylar commented May 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

brian-oneill commented May 8, 2026

Uh oh!

cbegeman commented May 8, 2026

Uh oh!

xylar commented May 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

hyungyukang left a comment

Choose a reason for hiding this comment

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

brian-oneill commented May 5, 2026 •

edited

Loading

Polaris `omega_pr` suite

xylar commented May 6, 2026 •

edited

Loading

Polaris `omega_pr` suite

cbegeman May 6, 2026 •

edited

Loading

xylar commented May 8, 2026 •

edited

Loading

xylar commented May 8, 2026 •

edited

Loading

xylar commented May 8, 2026 •

edited

Loading