Skip to content

Add velocity-component reconstruction for I/O - #525

Draft
xylar wants to merge 10 commits into
E3SM-Project:developfrom
xylar:omega/add-velocity-components
Draft

Add velocity-component reconstruction for I/O#525
xylar wants to merge 10 commits into
E3SM-Project:developfrom
xylar:omega/add-velocity-components

Conversation

@xylar

@xylar xylar commented Aug 25, 2026 •
edited
Loading

Copy link
Copy Markdown

Add VelocityZonalCell and VelocityMeridionalCell, the velocity components reconstructed at cell centers from the edge-normal velocity, as auxiliary fields that any IO stream can request.

Neither field is in any stream by default. They are computed once per time step, and only when a stream contains one of them. Stencil edges below the bathymetry count as zero velocity, as in MPAS-Ocean. On planar meshes the fields hold the Cartesian x and y components instead. This puts the VectorReconOnCell operator from #480 to use.

Other changes:

  • VectorReconOnCell can reconstruct one layer of a multi-layer field.
  • Vector reconstruction is available whenever the mesh file supplies the stencil and weights, not only on spheres. Planar meshes get this data from Fixes to planar vector reconstruction polaris#721, and the CTest planar mesh now has it.
  • SfcCoupling::init aborts if the mesh lacks the reconstruction data, which the coupler export needs.

Checklist

  • Documentation:

  • Linting

  • Building

    • CMake build does not produce any new warnings from changes in this PR

  • Testing

    aurora, oneapi-ifx, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    chrysalis, oneapi-ifx, openmpi

    • CTests Pass
    • Polaris omega_pr Pass

    frontier, craygnu, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    frontier, craygnu-mphipcc, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    pm-cpu, gnu, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    pm-gpu, gnugpu, mpich

    • CTests Pass
    • Polaris omega_pr Pass

  • Provide relevant details in a comment to the PR titled Testing with the following:

    • Which machines CTest unit tests
      have been run on and indicate that are all passing.
    • The Polaris omega_pr test suite
      has passed, using the Polaris e3sm_submodules/Omega baseline
    • Document machine(s), compiler(s), and the build path(s) used for -p for both the baseline (Polaris e3sm_submodules/Omega) and the PR build
    • Indicate "All tests passed" or document failing tests
    • Document testing used to verify the changes including any tests that are added/modified/impacted.

  • Performance related PRs: Please include a relevant PACE experiment link documenting performance before and after.

  • New tests:

    • CTest unit tests for new features have been added per the approved design.
    • Polaris tests for new features have been added per the approved design (and included in a test suite)

  • Stealth Features

    • If any stealth features are included in the PR, please confirm that they have been documented.

Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly.

@xylar

xylar commented Aug 25, 2026

Copy link
Copy Markdown
Author

This work is motivated by the need for reconstructed velocity components for analysis work in Polaris before the Sept. 15th deadline.

@xylar xylar changed the title Omega/add velocity components Add velocity-component reconstruction for I/O Aug 25, 2026
@xylar xylar mentioned this pull request Aug 25, 2026
Merged
2 tasks done
@xylar

xylar commented Aug 25, 2026

Copy link
Copy Markdown
Author

Needs to have the velocity reconstruction coefficients added to planar meshes for CTests, following the fix from E3SM-Project/polaris#721.

xylar added a commit to xylar/polaris that referenced this pull request Aug 25, 2026
@xylar
Use planar CTest mesh with reconstruction weights
11a8ace 
PlanarPeriodic48x48.omega_vars.260518.nc has no vector-reconstruction
fields, so Omega's planar CTests cannot exercise the reconstruction.
The two spherical CTest meshes gained them in the 260807 update, but
the planar one was left behind because the weights it would have gotten
were wrong.

Point at PlanarPeriodic48x48.omega_vars.260825.nc, which is the same
mesh with NEdgesReconOnCell, ReconStencilCell and ReconWeightsCell
appended, computed after the two planar fixes earlier on this branch.
No pre-existing variable changed. Omega will pick up the same datestamp
in E3SM-Project/Omega#525.
@xylar

xylar commented Aug 25, 2026

Copy link
Copy Markdown
Author

Early CTesting

CTest unit tests:

  • Machine: chrysalis
  • Compiler: intel
  • Build type: Release
  • Result: All tests passed
  • Log: /gpfs/fs1/home/ac.xylar/e3sm_work/polaris/fix-planar-velocity-reconstruction/build_omega/build_chrysalis_intel/ctests.log

@xylar
xylar requested a review from andrewdnolan August 25, 2026 21:41
xylar and others added 8 commits September 24, 2026 09:43
@xylar @claude
Add multi-layer support to VectorReconOnCell
104d2f9 
Add an overload that reconstructs one layer of a multi-layer field. The
reconstruction is independent in each layer, so this is the existing
single-layer form with a vertical index added.

Factor the Cartesian to local geographic rotation into a shared helper
so the spherical and planar paths stay identical between the two
overloads.

Extend the operator unit test to exercise the multi-layer form. Each
layer holds the exact field scaled by a layer-dependent factor so that
a reconstruction which mixed layers would show up as an error. The
error measures are normalized by the exact field, so the same expected
errors apply as for the single-layer case.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar @claude
Gate vector reconstruction on mesh contents, not geometry
cfbd29b 
The reconstruction stencil and weights are precomputed as a mesh
preprocessing step rather than by Omega, so whether they are available
is a property of the mesh file, not of whether the mesh is spherical.
Detect their presence by attempting the read of NEdgesReconOnCell and
treating a failure as "not present", and carry the answer as
Decomp::HasVectorRecon and HorzMesh::HasVectorRecon.

This lets planar meshes carry reconstruction data while leaving planar
meshes that predate it readable, and makes VectorReconOnCell abort with
a message naming the three arrays it needs rather than refusing every
planar mesh. The planar branch already present in the operator becomes
reachable as a result.

Decomp no longer parses the OnSphere attribute, which existed only to
gate these reads; HorzMesh still determines it from the full mesh
stream as before.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar @claude
Add VelocityZonalCell and VelocityMeridionalCell fields
cccbe35 
Add a VelocityReconAuxVars auxiliary variable class holding the zonal
and meridional velocity components reconstructed at cell centers from
the edge-normal velocity, and register both as fields in the
AuxiliaryState field group so they can be requested by any IO stream.
On a planar mesh the reconstructed vector already lies in the plane of
the mesh, so the two components are the Cartesian x and y components.

The reconstruction is written out here rather than calling
VectorReconOnCell so that layers outside MinLayerCell/MaxLayerCell keep
the fill value that Field::attachData wrote, which the operator, as a
reference implementation, knows nothing about.

Add AuxiliaryState::computeVelocityRecon to drive it. These are
diagnostics that no tendency reads, so the computation is deliberately
kept out of computeAll, which runs once per time stepper stage, and is
over owned cells only. It aborts if the mesh supplied no reconstruction
data. Nothing calls it yet.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar @claude
Compute the velocity components once per step, on demand
e3feaac 
Add IOStream::isFieldRequested, which reports whether a Field is in the
contents of any defined stream, so a diagnostic that nothing will read
or write can be skipped. Since validate replaces a group name in a
stream's contents by its member Fields, this is only meaningful after
the streams are validated; an unvalidated stream conservatively answers
yes so an optional computation is done rather than skipped.

Use it to gate AuxiliaryState::computeVelocityRecon, resolving the
answer on the first call and caching it, since the streams are
validated only after all Fields are defined and so well after the
auxiliary state is constructed.

Call it once per time step from both ocnRun loops, on the state the
step just produced (time level 0), which is the state that
updateTimeLevels has attached to the NormalVelocity field and so the
one that any stream writing this step will write.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar @claude
Check for vector reconstruction data in SfcCoupling::init
4caf3f9 
The surface velocity and SSH gradient exported to the coupler are
reconstructed at cell centers from edge-normal quantities, so a coupled
run needs a mesh that supplies the reconstruction stencil and weights.
Now that HorzMesh only reads these when the mesh file has them, check
HasVectorRecon at initialization so a coupled run on a mesh without them
fails there rather than at the first export.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
@xylar @claude
Document the reconstructed velocity components
08dad9c 
Add VelocityZonalCell and VelocityMeridionalCell to the auxiliary
variable list in the user guide, with a section on the three ways they
differ from the other auxiliary variables: they are diagnostic and so
computed once per step, they are only computed when a stream asks for
them, and they need a mesh file that supplies the reconstruction stencil
and weights.

Document computeVelocityRecon in the auxiliary state developer guide,
and IOStream::isFieldRequested in the IOStreams developer guide.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar @claude
Point the developer guide at the new planar test mesh
a1e03fc 
PlanarPeriodic48x48.omega_vars.260825.nc is the first planar mesh to
carry the vector reconstruction stencil and weights, so the CTests that
use a planar mesh need it rather than the earlier datestamps.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar @claude
Test vector reconstruction on planar meshes
cf86cf3 
Drop the early return that made testVectorRecon a no-op for planar
meshes, now that the planar test mesh carries the reconstruction
stencil and weights. On a plane the reconstructed components are the
Cartesian x and y ones rather than zonal and meridional, which the
magnitude comparison the test makes is indifferent to.

The expected error norms for the planar setup were derived offline from
the mesh file, using a calculation that reproduces the spherical
expected values already in this file to within 5e-14 relative. They
have since been confirmed by a CTest run on chrysalis with intel.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar
xylar force-pushed the omega/add-velocity-components branch from b9da99c to cf86cf3 Compare September 24, 2026 14:55
@xylar

xylar commented Sep 24, 2026 •
edited
Loading

Copy link
Copy Markdown
Author

Testing

On chrysalis with intel and openmpi:

  • CTests at 9bea205: all 50 pass (with the Omega#572 enable_language(Fortran) workaround). The new velocity reconstruction check in AUXVARS_PLANE_TEST and AUXVARS_SPHERE_TEST fails without 94306e2.
  • omega_pr at 94306e2 (the last source change): all pass, bit-for-bit against the Polaris e3sm_submodules/Omega baseline. Baseline build: /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260921/build_omega_develop_5dfd2ae971, PR build: /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260924/velrecon525_fix/ocean_build
  • With VelocityZonalCell and VelocityMeridionalCell added to History in cosine_bell, baroclinic_channel and QU240 realistic_global, they match an independent Python reconstruction to round-off. cosine_bell is within 7e-4 of the analytic solid-body velocity and identical on 12 and 24 ranks.

Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly.

xylar and others added 2 commits September 24, 2026 10:46
@xylar @claude
Skip inactive stencil edges in the velocity reconstruction
94306e2 
The reconstruction stencil of a cell reaches beyond its own edges, so
near bathymetry it includes edges that are not active in a layer where
the cell is. Their normal velocity is zero or, below both neighboring
cells, the fill value, which leaked into the reconstructed components.
On the QU240 realistic_global mesh this gave values near 1e35 in 7% of
the active cell layers.

Treat edges outside [MinLayerEdgeBot, MaxLayerEdgeTop] as having zero
normal velocity, as MPAS-Ocean's zeroed normalVelocity does there.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
@xylar @claude
Test the velocity reconstruction next to bathymetry
9bea205 
Add a VelocityReconAuxVars check to AuxiliaryVarsTest in which every
fourth column (by global cell ID, so halos agree) is half as deep. The
velocity is masked as the model masks it, with the fill value below both
cells of an edge and zero below one, and the result is compared in each
cell's active layers against VectorReconOnCell applied to the same
velocity with every inactive edge set to zero.

Without the previous commit this fails on both meshes, with differences
near 1e36 from the fill value. It runs last since it changes the
vertical layer ranges the other checks assume are full depth.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
@xylar

xylar commented Sep 25, 2026 •
edited
Loading

Copy link
Copy Markdown
Author

Testing

On pm-gpu with gnugpu and mpich:

  • CTests at 9bea205: all 50 pass. The Omega#572 enable_language(Fortran) workaround is not needed for this compiler.
  • omega_pr at 9bea205: all 25 tasks pass, bit-for-bit against the Polaris e3sm_submodules/Omega baseline across 141 comparisons. Baseline build: /pscratch/sd/x/xylar/polaris_1.1/pm-gpu/test_20260924/velrecon525/baseline/ocean_build, PR build: /pscratch/sd/x/xylar/polaris_1.1/pm-gpu/test_20260924/velrecon525/pr/ocean_build
  • With VelocityZonalCell and VelocityMeridionalCell added to the History contents in cosine_bell and baroclinic_channel (realistic_global already writes the AuxiliaryState group), they match an independent Python reconstruction built from the mesh stencil and weights:
task ranks max relative difference
icos/cosine_bell/decomp 12, 24 5.5e-16
baroclinic_channel/10km/decomp 4, 8 2.2e-16
QU.240km/analysis_members_test 1 2.8e-16

Both fields are bit-for-bit between the two rank counts in each decomp task, and cosine_bell is within 7.0e-4 of the analytic solid-body velocity. QU240 leaves 16% of cell-layers inactive, with columns from 2 to 60 layers, so the stencil's inactive-edge skip is exercised there.

Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly.

@xylar xylar added enhancement New feature or request Omega labels Oct 21, 2026
@vanroekel vanroekel added this to the Feature Freeze milestone Oct 24, 2026

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@andrewdnolan andrewdnolan Awaiting requested review from andrewdnolan

Assignees

No one assigned

Labels

enhancement New feature or request Omega

Projects

None yet

Milestone

Feature Freeze

Development

Successfully merging this pull request may close these issues.

2 participants

@xylar @vanroekel