Skip to content

Add centered pressure gradient implimentation#8

Open
sbrus89 wants to merge 159 commits intoomega/developfrom
omega/pgrad
Open

Add centered pressure gradient implimentation#8
sbrus89 wants to merge 159 commits intoomega/developfrom
omega/pgrad

Conversation

@sbrus89
Copy link
Copy Markdown
Owner

@sbrus89 sbrus89 commented Oct 17, 2025

No description provided.

Copy link
Copy Markdown
Collaborator

@xylar xylar left a comment

Choose a reason for hiding this comment

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

@sbrus89, a few comments. The main ones in terms of needing to discuss them relate to how the product of $\alpha p$ should be computed at layer interfaces and edges. I think it probably needs to be computed based on 4-point averages from cell centers and layer midpoints but we should make sure we're on the same page.

Comment thread components/omega/src/ocn/Eos.cpp Outdated
parallelFor(
"compute-interface-values", {NCellsAll, NVertLayers + 1},
KOKKOS_LAMBDA(I4 ICell, I4 KInterface) {
if (KInterface == 0) {
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

As we discussed, we need hierarchical parallelism here and this should be MinLayerCell or whatever the proper variable is, rather than 0.

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

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

Yes, I'll take care of swapping out for hierachical parallelism.

Comment thread components/omega/src/ocn/Eos.cpp Outdated
// AbsSalinityInterface(ICell, KInterface) =
// AbsSalinity(ICell, KInterface);
SpecVolInterface(ICell, KInterface) = SpecVol(ICell, KInterface);
} else if (KInterface == NVertLayers + 1) {
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Similarly, this should have hierarchical parallelism and use MaxLayerCell or similar.

Comment thread components/omega/src/ocn/Eos.cpp Outdated
Comment thread components/omega/src/ocn/PGrad.h Outdated
const I4 ICell1 = CellsOnEdge(IEdge, 1);
const Real InvDcEdge = 1.0_Real / DcEdge(IEdge);

for (int KVec = 0; KVec < VecLength; ++KVec) {
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This needs to be over something like MinLayerEdgeBot to MaxLayerEdgeTop with hierarchical parallelism.

Comment thread components/omega/src/ocn/PGrad.h Outdated
Comment on lines +53 to +54
(PressureInterface(ICell1, K) * SpecVolInterface(ICell1, K) +
PressureInterface(ICell0, K) * SpecVolInterface(ICell0, K));
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Let's chat about it but I think we want the PAlpha product to be averaged both to cell centers and to layer interfaces.

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

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

Check out my latest commit and let me know what you think. The new averaging appeared to lower errors in my hydrostatic test.

sbrus89 pushed a commit that referenced this pull request Dec 3, 2025
…t#7866)

This is changing code that is widely used and has been unchanged for a
very long time, so I'm adding a lot of reviewers. Feel free to add
more.

I stumbled across this code when trying to debug a weird eamxx failure:

I observed the following symtoms:

If you look at components/elm/src/main/reweightMod.F90, you'll see the line SHR_ASSERT(bounds%level == BOUNDS_LEVEL_CLUMP, errMsg(__FILE__, __LINE__)) which causes this error:

1: free(): invalid pointer
...
0: #8  0x4c7e3ab in __shr_log_mod_MOD_shr_log_errmsg
0:      at /pscratch/sd/a/acmetest/E3SM/share/util/shr_log_mod.F90:78
0: E3SM-Project#9  0x914149 in __reweightmod_MOD_reweight_wrapup
0:      at /pscratch/sd/a/acmetest/E3SM/components/elm/src/main/reweightMod.F90:48
0: E3SM-Project#10  0x17bda7d in __dynsubgriddrivermod_MOD_dynsubgrid_wrapup_weight_changes
0:      at /pscratch/sd/a/acmetest/E3SM/components/elm/src/dyn_subgrid/dynSubgridDriverMod.F90:403
0: E3SM-Project#11  0x17bf256 in __dynsubgriddrivermod_MOD_dynsubgrid_init._omp_fn.0
0:      at /pscratch/sd/a/acmetest/E3SM/components/elm/src/dyn_subgrid/dynSubgridDriverMod.F90:150

I thought this was because the SHR_ASSERT was failing, but that's not
the case (both bounds%level and BOUNDS_LEVEL_CLUMP are always 2). If I
just comment this SHR_ASSERT out entirely, the test PASSES!

The lowest line of code that triggers the free(): invalid pointer error is:
shr_log_errMsg = 'ERROR in '//trim(file)//' at line '//toString(line)

What I think is happening is that there is some memory corruption that
is causing the freeing of the allocation done by toString to fail. So,
this PR kind of sweeps that issue under the rug, but I do think it's
better not to do dynamic allocations if you don't have to. The test
PASSes with the changes in this PR.

[BFB]
@xylar
Copy link
Copy Markdown
Collaborator

xylar commented Dec 4, 2025

@sbrus89, your last few commits look reasonable to me. I would remove the commented-out loops after the switch to hierarchical parallelism eventually but it's okay to keep them for now while you make sure things are working as expected.

It's getting cumbersome to review here because your local omega/develop is out of date, so let's just got ahead and make the PR to Omega.

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

// compute specific volume from EOS
VCoord->computePressure(LayerThick, SurfacePressure);
DefEos->computeSpecVol(Temp, Salinity, PressureMid);
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This is probably not where we want to be computing the specific volume in the long run, right? We would want to compute it somewhere external in the time step and just make it available here? I guess that's what Hyun is working on?

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

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

Right, I realized I need to add this call in Tendencies as well.

Comment thread components/omega/src/ocn/PGrad.h Outdated
Comment on lines +49 to +50
Real Rho0 = 1026.0_Real;
Real Gravity = 9.80616_Real;
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This is not the right Gravity. We should use the value from GlobalConstants.h, which comes from here:
https://github.com/E3SM-Project/Omega/blob/develop/share/util_cxx/pcd_const.h#L28

Comment thread components/omega/src/ocn/PGrad.h Outdated
const I4 ICell0 = CellsOnEdge(IEdge, 0);
const I4 ICell1 = CellsOnEdge(IEdge, 1);
const Real InvDcEdge = 1.0_Real / DcEdge(IEdge);
Real Rho0 = 1026.0_Real;
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This isn't used and should be removed.

mwarusz and others added 8 commits February 20, 2026 11:32
Surface flux terms, body forces and pressure terms were all getting
multiplied by $\rho_0$, leading to incorrect dimensions and
inconsistency with the rest of the derivation.
…1-equations

Fixes to Omega v1 Governing Equations

This pull request updates the documentation in OmegaV1GoverningEqns.md to improve the consistency, clarity, and correctness of the mathematical notation used in the derivation of the governing equations. The changes primarily remove unnecessary factors of \rho_0 from several equations, clarify the definition and notation for surface forces, and correct integration limits and variable usage for consistency.

Mathematical and Notational Corrections:

* Removed unnecessary factors of \rho_0 from the pseudo-height mass, tracer, and momentum conservation equations, as well as from their layer-averaged and Reynolds-averaged forms, to ensure dimensionally consistent and simplified expressions.
* Clarified the definition of the surface force (traction vector) as \tau = \sigma \cdot {\bf n} and distinguished the horizontal component \tau_\perp in the momentum equation, improving the physical interpretation and accuracy of boundary force terms.

Consistency and Integration Limit Corrections:

* Corrected integration limits and variable usage in the layer and Reynolds-averaged tracer equations to ensure consistency with the rest of the document and the intended physical meaning.
* Updated several equations to use consistent notation for averages, fluctuations, and layer integrals, improving readability and reducing the risk of confusion.
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.

8 participants