Skip to content

Add MOC stream function analysis capability - #481

Merged
xylar merged 53 commits into
E3SM-Project:developfrom
brian-oneill:omega/analysis-moc-streamfunction
Sep 28, 2026
Merged

xylar merged 53 commits into
E3SM-Project:developfrom
brian-oneill:omega/analysis-moc-streamfunction

Conversation

@brian-oneill

@brian-oneill brian-oneill commented Jul 28, 2026 •

Copy link
Copy Markdown

Overview

This PR introduces a complete MOC (Meridional Overturning Circulation) analysis capability to Omega, enabling computation of the MOC streamfunction using two complementary methods:

  1. Latitude-binned regional MOC: Computes MOC as a function of latitude and depth for specified ocean regions
  2. Transect-based MOC: Computes MOC across specific transects as a function of depth

The implementation adds 8 new analysis operators, enhanced analysis infrastructure with regional mask support, a new MOC analysis group, and comprehensive configuration templates.

Key Features

New Analysis Operators (8 total):

  • BinaryMultiplyOp: Element-wise field multiplication with vertical expansion support
  • BinnedAccumulatorOp: Accumulates field values into spatial bins (core MOC operator)
  • CoordinateBinningOp: Assigns mesh entities to bins based on coordinate values
  • ExtractRegionOp: Applies regional masks to fields
  • PrefixSumOp: Cumulative summation (integration) along specified dimension
  • PseudoToGeometricOp: Converts pseudo-height quantities to geometric coordinates
  • ScalarMultiplyOp: Multiplies field by scalar constant for unit conversion
  • TransectAccumulatorOp: Accumulates transport across transect edges

Infrastructure Enhancements:

  • Field Class Regional Mask Support: Fields can carry spatial mask information through operator chains with automatic propagation
  • IOStream IOName Metadata: Provides user-friendly variable names in netCDF output while maintaining internal operator chain naming
  • Enhanced AnalysisGroup Base Class: New setOutputIOName method and operator-specific configuration support

MOC Analysis Group:

  • Bundled analysis group for computing MOC streamfunction
  • Configurable latitude binning (number of bins, lat range)
  • Regional MOC computation for named ocean regions
  • Transect-based MOC computation across specified transects
  • Configurable temporal output (reduction periods and snapshots)
  • IOStream integration for netCDF output

MOC Computation Pipeline

Latitude-binned Regional MOC chain:

  1. Assign cells to latitude bins (static, initialization)
  2. Convert vertical pseudo-velocity to geometric coordinates
  3. Compute vertical flux (velocity × area)
  4. Apply regional mask (optional)
  5. Accumulate into latitude bins
  6. Horizontal integration (south→north)
  7. Convert to Sverdrups

Transect-based MOC chain:

  1. Convert pseudo-thickness to geometric layer thickness
  2. Compute transport (thickness × velocity × edge width)
  3. Accumulate across transect edges
  4. Vertical integration (bottom→top)
  5. Convert to Sverdrups

Technical Implementation

Design Features:

  • Field-level regional mask storage and automatic propagation
  • IOName metadata system separating internal from user-facing names
  • Kokkos hierarchical parallelism for 2D/3D operators
  • SFINAE compile-time optimization to prevent invalid template instantiations
  • Operator parameters passed via Config objects through chain parsing

Output Format:

  • MOC streamfunction in Sverdrups (1 Sv = 10⁶ m³/s)
  • Dimensions: latitude bins × depth levels (regional), or depth only (transect)
  • NetCDF files with user-friendly variable names (e.g., "MOC_streamfunction_Global")
  • Configurable temporal averaging and snapshot output

Limitations

  • Regional masks not yet implemented (operator infrastructure ready, placeholder in config)
  • Transect masks not yet implemented (operator infrastructure ready, placeholder in config)

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.
  • 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)

@brian-oneill
brian-oneill requested review from cbegeman and xylar July 28, 2026 04:01
@xylar

xylar commented Jul 28, 2026

Copy link
Copy Markdown

Sorry, @brian-oneill, I didn't get to this today. I'll try again tomorrow.

Also, let me know what you need from me regarding both the dynamic streams here in Omega and the Polaris support.

Comment on lines +405 to +406
ReductionPeriod: [1Month]
SnapshotPeriod: [1Day]

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Can you help me understand what sets how often the MOC is computed? Is it the SnapshotPeriod? With the options above, ReductionPeriod would then be averaging daily instantaneous MOC values over 1 month?

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

My understanding is that they are independent of one another. In practice, we would have:

  ReductionPeriod: [1Month]
  SnapshotPeriod: []

since we want monthly averages and don't need snapshots.

It's hard for me to imagine very much analysis where we want both time averages and snapshots at the same time, in practice.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

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

Yes, ReductionPeriod and SnapshotPeriod produce outputs independently. The ReductionPeriod outputs currently accumulate every timestep so the MOC is computed each timestep for time averages, but this can be extended to allow for a courser sampling frequency pretty easily.

# Computes spatial reduction statistics (Mean, Min, Max, StdDev)
# for a set of ocean fields. Supports temporal reduction (time-averaged
# output over a window) and instantaneous snapshots (discrete sampling).
Fields: [NormalVelocity, PseudoThickness, Temperature, Salinity]

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Do we need to test that when, e.g., LayerThickness_BinaryMultiply(NormalVelocity), is present here that the MOC chain uses the available field or is this kind of thing covered by existing CTests?

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Well, I guess not present here because GlobalStats reduces spatially.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

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

During the initial parsing, the parser checks if a Field that would be output by an operator has already been registered, to prevent building a duplicate operator. There is a unit test that checks this behavior, but it could be more robust.

@xylar xylar left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

@brian-oneill, I'll do some testing but here are a few comments to keep the process moving.

This looks great! A lot of the pieces are in place and just a few tweaks would be helpful, I think. Plus a few things that might be for now or might be postponed until later.

Comment thread components/omega/configs/Default.yml Outdated
Comment on lines +405 to +406
ReductionPeriod: [1Month]
SnapshotPeriod: [1Day]

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

My understanding is that they are independent of one another. In practice, we would have:

  ReductionPeriod: [1Month]
  SnapshotPeriod: []

since we want monthly averages and don't need snapshots.

It's hard for me to imagine very much analysis where we want both time averages and snapshots at the same time, in practice.

Comment thread components/omega/src/analysis/analysisGroups/MOC.cpp
Comment thread components/omega/src/analysis/operators/PseudoToGeometricOp.h
Comment thread components/omega/src/analysis/operators/PseudoToGeometricOp.h

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Could ScalarMultiply be generalized to take a model config option or known constant as its input, not just a hard-coded number? This would seem much more useful and general.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

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

For the use case here, this Op gets passed 1e-6 through a "config" that is programmatically defined in the parser. So 1e-6 is hard-coded into the MOC chain construction, but the Op itself is designed to take a configurable value to be compatible with the future composable framework.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Sounds good. It's RhoSw in particular that I wanted to know about.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

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

Gotcha, that's possible but requires a little more work. Allowing users to request a defined constant from the config file would require defining a map between string labels and the variable names in GlobalConstants.h

std::map<std::string,Real> Constants = {
{"RhoSw", RhoSw},
{"Gravity", Gravity},
{"Pi", Pi},
...
};

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I think we should try to figure out how to automate that somehow but, yes, that's what I was anticipating. Nothing that needs to be in this PR.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

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

Yeah, automated would be best. Here's Claude's suggestion:

Yes, the X-macro pattern is exactly right. The idea: define a single list of constants once using a macro, then expand it in two ways — once to declare the constexpr values, and once to build the lookup.

In GlobalConstants.h, replace the individual constexpr declarations for the physical constants with:

// X-macro list: X(Name, Value)
#define OMEGA_PHYSICAL_CONSTANTS(X)          \
    X(RhoSw,   pcd::seawater_density_reference)           \
    X(RhoFw,   pcd::pure_water_density_reference)         \
    X(RhoAir,  pcd::dry_air_density_at_standard_temperature_and_pressure) \
    X(Gravity, pcd::standard_acceleration_of_gravity)     \
    X(RhoIce,  pcd::sea_ice_density_reference)            \
    /* ... all others ... */

// Expand to constexpr declarations (same as before)
#define DECLARE_CONST(Name, Value) constexpr Real Name = (Value);
OMEGA_PHYSICAL_CONSTANTS(DECLARE_CONST)
#undef DECLARE_CONST

Then the lookup function writes itself:

inline std::optional<Real> getConstantByName(const std::string &Name) {
#define MATCH_CONST(CName, Value) if (Name == #CName) return CName;
    OMEGA_PHYSICAL_CONSTANTS(MATCH_CONST)
#undef MATCH_CONST
    return std::nullopt;
}

#CName stringifies the identifier automatically, so the name in the config file ("RhoSw") matches the C++ variable name without any manual duplication.

Pros: Zero maintenance burden — adding a new constant to the list automatically makes it available by name to ScalarMultiplyOp and any future operator that calls getConstantByName.

Cons: Requires refactoring the existing constexpr declarations in GlobalConstants.h into the macro list format. The math-derived ones (TwoPi, SDay, etc.) are trickier since they depend on other constants — those would need to stay as regular constexpr or be added after the macro expansion.

The practical approach: put only the "leaf" physical constants (densities, heat capacities, etc.) in the X-macro list, and keep derived/compound ones (TwoPi, TkFrzSw, SDay) as regular constexpr declarations below. You could optionally add those to a second X-macro list if you want them accessible by name too.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Not bad, not bad!

@xylar xylar left a comment •

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I successfully ran a 5-day test with the MOC on in:
/lcrc/group/e3sm/ac.xylar/polaris_1.0/chrysalis/test_20260804/ec30to60-global-moc/ocean/spherical/realistic_global/EC30to60E2r2/analysis_members_test

I used:

    MOC:
      Enable: true
      # Meridional Overturning Circulation (MOC) streamfunction analysis group
      # Computes MOC as a function of latitude and depth for regions,
      # and as a function of depth for transects
      NumBins: 180              # Number of latitude bins (default: 180, ~1 degree)
      MinLat: -90.0             # Minimum latitude in degrees (default: -90.0)
      MaxLat: 90.0              # Maximum latitude in degrees (default: 90.0)
      Regions: [Global]         # List of region names for regional MOC
                                # NOTE: Region masks not yet implemented
      Transects: []             # List of transect names for transect-based MOC
                                # NOTE: Transect masks not yet implemented
      ReductionPeriod: [1day]       # Temporal reduction periods
      SnapshotPeriod: []  # Instantaneous output periods
      Filename: moc.$Y
      Stream:
        FileFreq: 1
        FileFreqUnits: days

So daily averaging rather than monthly for efficiency.

Here's an example plot:
Image

However, the latitude bins and depth are missing from the output file:

$ ncdump -h moc_1dayTimeStats.0001 
netcdf moc_1dayTimeStats {
dimensions:
	MaxCellsOnEdge = 2 ;
	MaxEdges = 7 ;
	MaxEdges2 = 14 ;
	NCells = 236853 ;
	NEdges = 719506 ;
	NTracers = 2 ;
	NVertLayers = 60 ;
	NVertLayersP1 = 61 ;
	NVertices = 482371 ;
	NumBinsLatCell_BinIndex = 180 ;
	Scalar = 1 ;
	VertexDegree = 3 ;
	time = UNLIMITED ; // (5 currently)
variables:
	double MOC_streamfunction_Global_TimeMean1day(time, NumBinsLatCell_BinIndex, NVertLayersP1) ;
		MOC_streamfunction_Global_TimeMean1day:Description = "Time average of VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)" ;
		MOC_streamfunction_Global_TimeMean1day:Name = "VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)_TimeMean1day" ;
		MOC_streamfunction_Global_TimeMean1day:StdName = "" ;
		MOC_streamfunction_Global_TimeMean1day:Units = "" ;
		MOC_streamfunction_Global_TimeMean1day:ValidMax = 1.79769313486232e+308 ;
		MOC_streamfunction_Global_TimeMean1day:ValidMin = -1.79769313486232e+308 ;
		MOC_streamfunction_Global_TimeMean1day:_FillValue = 9.96920996838687e+36 ;
		MOC_streamfunction_Global_TimeMean1day:long_name = "Time average of VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)" ;
		MOC_streamfunction_Global_TimeMean1day:name = "VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)_TimeMean1day" ;
		MOC_streamfunction_Global_TimeMean1day:standard_name = "" ;
		MOC_streamfunction_Global_TimeMean1day:units = "" ;
		MOC_streamfunction_Global_TimeMean1day:valid_max = 1.79769313486232e+308 ;
		MOC_streamfunction_Global_TimeMean1day:valid_min = -1.79769313486232e+308 ;
	double time(time) ;
		time:Description = "time" ;
		time:Name = "time" ;
		time:StdName = "time" ;
		time:Units = "seconds since 0001-01-01 00:00:00" ;
		time:ValidMax = 1.e+20 ;
		time:ValidMin = 0. ;
		time:_FillValue = 9.96920996838687e+36 ;
		time:calendar = "noleap" ;
		time:long_name = "time" ;
		time:name = "time" ;
		time:standard_name = "time" ;
		time:units = "seconds since 0001-01-01 00:00:00" ;
		time:valid_max = 1.e+20 ;
		time:valid_min = 0. ;

// global attributes:
		:SimulationTime = "0001-01-06_00:00:00" ;
		:SimulationTime0 = "0001-01-02_00:00:00" ;
		:SimulationTime1 = "0001-01-03_00:00:00" ;
		:SimulationTime2 = "0001-01-04_00:00:00" ;
		:SimulationTime3 = "0001-01-05_00:00:00" ;
		:SimulationTime4 = "0001-01-06_00:00:00" ;
}

Also the file is missing a .nc extension.

Happy to rerun once this is fixed.

@brian-oneill

brian-oneill commented Aug 4, 2026 •

Copy link
Copy Markdown
Author

Is refBottomDepth the appropriate depth to go with the output? Because that's not currently read in to Omega...

I've noticed the history and hifreq files that get created when running the omega ctests also no longer get the .nc suffix added. Looks like 1f92657 removed the lines that explicitly added the suffix to filenames with time templates. A comment says PIO should handle this, but that doesn't appear to be the case. Adding .nc to the end of the Filename option in the config does work.

@xylar

xylar commented Aug 4, 2026

Copy link
Copy Markdown

refBottomDepth is an MPAS-Ocean concept and is not going to be available in general for Omega. So, no, it is not appropriate.

The problem with the current layer-wise approach is it implicitly assumes pure z-level layers. The best way to get depths given this approach is to get the area-weighted average of zInterface and use that as the vertical coordinate. For now, that's fine. With ice-shelf cavities, sigma coordinates, etc. in the future, we'll need to do vertical binning or interpolation to a z-level or density-level grid.

@xylar

xylar commented Aug 4, 2026 •

Copy link
Copy Markdown

I'm fine with whatever fix to the .nc extension but we do need them to be put back somehow. I can understand the desire to support arbitrary format and not explicitly require NetCDF, though.

# List of field names to compute statistics for
SpatialStats: [Max, Min, Mean, StdDev]
# Spatial statistics to compute (one per field)
ReductionPeriod: [1Day, 1Month]

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Suggested change
ReductionPeriod: [1Day, 1Month]
ReductionPeriod: []

Should we remove reductions altogether from this analysis member to save compute time?

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I think we want 1Month reduction for the climatology plot, don't we?

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I guess it depends on what configuration Default.yml is aimed at -- typical standalone testing or a longer production run.

# Spatial statistics to compute (one per field)
ReductionPeriod: [1Day, 1Month]
# Temporal reduction periods (time-averaged stats)
SnapshotPeriod: [6Hours]

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

MPAS-O is equivalent to 1Day but I like increasing the freq to 6 hours

@cbegeman

Copy link
Copy Markdown

However, the latitude bins and depth are missing from the output file:

@brian-oneill Is this ready for testing with the fix to this comment?

@brian-oneill

Copy link
Copy Markdown
Author

However, the latitude bins and depth are missing from the output file:

@brian-oneill Is this ready for testing with the fix to this comment?

@cbegeman I added the latitude bins. The depths will require a bit more thought and work, and I won’t have time to complete that piece until I’m back at the end of the week.

@cbegeman

cbegeman commented Sep 1, 2026 •

Copy link
Copy Markdown

Polaris omega_pr suite

  • Baseline workdir: //lcrc/group/e3sm/ac.cbegeman/scratch/MPAS-Ocean-output/chrys/polaris-main-omega-submodule-20260901/
  • Baseline build: /lcrc/group/e3sm/ac.cbegeman/scratch/MPAS-Ocean-output/chrys/polaris-main-omega-submodule-20260901/build
  • PR build: /lcrc/group/e3sm/ac.cbegeman/scratch/MPAS-Ocean-output/chrys/polaris-main-omega-moc-20260901/build
  • PR workdir: /lcrc/group/e3sm/ac.cbegeman/scratch/MPAS-Ocean-output/chrys/polaris-main-omega-moc-20260901
  • Machine: chrysalis
  • Partition: debug
  • Compiler: intel
  • Build type: Release
  • Log: /lcrc/group/e3sm/ac.cbegeman/scratch/MPAS-Ocean-output/chrys/polaris-main-omega-moc-20260901/polaris_omega_pr.o1279292
  • Result: All tests passed

@cbegeman

cbegeman commented Sep 3, 2026

Copy link
Copy Markdown

@brian-oneill I successfully ran a 10-day QU240 realistic_global test with this feature enabled and daily time reduction. I inspected the moc file's attributes and variables but didn't plot them, given that the latest changes would not have affected the MOC computation itself. GeomZInterface_HorzMean looks reasonable and only varies slightly over my 10-day run as expected.

@cbegeman cbegeman left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Approving on the basis of visual inspection, my testing on chrysalis-intel, and @xylar testing and viz. Great work, @brian-oneill !

@xylar

xylar commented Sep 7, 2026

Copy link
Copy Markdown

@brian-oneill, I'm reviewing but it's turning up a need to rebase this branch onto develop to get CIME updates. Could you do the rebase when you have a chance. It seems to be clean.

@xylar xylar left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

@brian-oneill, this is great and I think we're almost there!

How this was reviewed

I reviewed this by running the omega_analysis suite from the Polaris branch in E3SM-Project/polaris#743. That PR adds an ocean/analysis/moc step to Polaris' Omega analysis suite: it reads this group's output, averages the reductions over a range of years weighted by the length of each period, and plots the streamfunction against latitude and the interface elevations the group provides. Everything below came out of getting that to work.

The simulation. QU240, one simulated year from 0000-12-01, with

MOC:
  Enable: true
  NumBins: 60
  MinLat: -90.0
  MaxLat: 90.0
  Regions: [Global]
  Transects: []
  ReductionPeriod: [1Month]
  SnapshotPeriod: []
  Filename: moc.$Y-$M.nc
  Stream: {FileFreq: 1, FileFreqUnits: months}

built from this branch rebased onto omega/develop, on chrysalis with intel (oneAPI 2025.2).

The analysis. Polaris' omega_analysis suite, all five tasks passed.

The plot, and the rest of the suite's products, are on the LCRC portal — the MOC is the moc_global_0001-0001 gallery:

https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/xasaydavis/omega_analysis_moc_20260907/

Work directories: the run at /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260907/qu240-moc-analysis-mockup, the analysis at .../test_20260907/qu240-moc-analysis.

What the previous review asked for, and what is now there

The previous review reported three things: latitude bins missing, depths missing, and no .nc extension. All three are addressed.

double MOC_streamfunction_Global_TimeMean1Month(time, NumBinsLatCell_BinIndex, NVertLayersP1) ;
double MOCLatBinBoundaries(time, NMocLatBinBoundaries) ;
double GeomZInterface_HorzMean(time, NVertLayersP1) ;
double time(time) ;

NumBinsLatCell_BinIndex = 60, NMocLatBinBoundaries = 61, NVertLayersP1 = 61. The file is moc_1MonthTimeStats.0001-01.nc, so the suffix is there too.

What works

Lots is working exactly as needed here!

The streamfunction closes exactly. In every month of the run, psi is exactly zero at the surface interface, exactly zero at the bottom interface, and exactly zero in the southernmost bin. The residual at the northernmost bin — the one place a global streamfunction has anything left to accumulate — is small and shrinking across the run (0.255, 0.138, 0.115 Sv over the first three months). The bottom-interface zero is the check the earlier cross-validation script used, and it holds to round-off.

Reductions are genuinely per-period. Consecutive monthly means differ substantially (up to 46 Sv between months 2 and 3), so the temporal accumulator is being reset between periods rather than producing a running average.

Reductions are stamped at the end of their period. The time values are 31, 62, 90 and 121 days since 0000-12-01, which are the ends of December, January, February and March. That makes the period each reduction covers recoverable from the stamps.

The file names are reconstructable from the config alone, as <prefix>_<period><TimeStats|Instants><template>, and the variable names as MOC_streamfunction_<Region>_TimeMean<period>. Polaris builds both from the simulation's own YAML without any of it being restated, which is what makes the analysis configuration-free.

Findings

1. MOCLatBinBoundaries should very likely not have a time dimension

It is static: computed once in the MOC constructor from NumBins, MinLat and MaxLat, and never updated. But it is attached to the output streams with addField() like any other field, so every reduction in every file carries a copy of the same 61 numbers.

The cost in bytes is nothing. The cost to a consumer is not: a reader that takes the variable as written gets NumBins + 1 boundaries per output time, and over a twelve-month range that flattens to twelve times as many latitudes as there are bins. Polaris now takes the first record and checks the rest agree, but that is a special case that exists only because the field claims to vary and does not.

GeomZInterface_HorzMean having a time dimension is correct — the interfaces move — so this is specifically about the bin boundaries.

2. The streamfunction still carries no units

MOC_streamfunction_Global_TimeMean1Month:Units = "" ;
MOC_streamfunction_Global_TimeMean1Month:units = "" ;

The chain ends in ScalarMultiply(1.0e-6), so the values are in Sverdrups, but nothing in the file says so and a consumer has to know it out of band. This was raised in the previous review and is unchanged. Both spellings need Sv: a PR addressing #529 will remove the capitalized Units in favor of the CF-compliant units, but until it lands both have to be correct.

3. GeomZInterface_HorzMean fills in only one of its two units attributes

GeomZInterface_HorzMean:Units = "" ;
GeomZInterface_HorzMean:units = "m" ;

The duplicated attributes are already tracked in #529; a fix would drop the capitalized set in favor of the CF-compliant ones. Until that lands both have to be right, and here only the lower-case one is, so a reader picking Units gets nothing. The same doubling appears on Name/name, StdName/standard_name and ValidMin/valid_min.

4. The written bin boundaries are not the edges the operator bins on

MOC.cpp writes the boundaries without the margin that CoordinateBinningOp applies before it bins:

// MOC.cpp
const Real BinWidthDeg = (MaxLat - MinLat) / static_cast<Real>(NumBins);
BinBoundHost(I) = MinLat + I * BinWidthDeg;
// CoordinateBinningOp.h
const Real Margin = (MaxBin - MinBin) * 1.0e-6;
MinBin -= Margin;
MaxBin += Margin;
BinWidth = (MaxBin - MinBin) / NumBins;

For the defaults the difference is about 1.8e-4 degrees, far below anything that matters for a plot. Even so, it's better write out what the code really uses. Deriving the written boundaries from the operator's own MinBin/BinWidth would make them agree by construction. Alternatively, we should build the margin into the input latitude range somehow.

5. Reductions carry no time bounds

The output has a CF time with units and calendar but no time_bnds, and no cell_methods saying the values are a time mean. A consumer averaging reductions over a range has to weight each by the length of its period — the streamfunction is linear in the velocity, so this is not optional, and monthly periods differ in length by up to a tenth. With bounds that is exact; without them Polaris infers each period from the gap to the next stamp, which is exact for every period except the first, whose length it has to assume.

This is not specific to the MOC group; it applies to any temporal reduction, and it is the same gap that makes ncclimo unusable on Omega history files without a pre-processing step.

6. The long_name is the operator chain

MOC_streamfunction_Global_TimeMean1Month:long_name =
  "Time average of VerticalPseudoVelocity_PseudoToGeometric_BinaryMultiply(AreaCell)_BinnedAccumulator(LatCell_BinIndex)_PrefixSum_ScalarMultiply(1.0e-6)" ;

Excellent provenance and exactly what Description ( renamed to comment when #529 gets addressed) should say. As a long_name it is what ends up on a colorbar or in a variable listing, where "global meridional overturning streamfunction" would serve better.

7. Noted along the way

Running the MOC group with a monthly ReductionPeriod aborts unless RestartWrite uses a Months or Years interval, because AnalysisGroup validates the restart interval without checking that RestartWrite has a periodic alarm. Filed separately as #539.

@brian-oneill
brian-oneill force-pushed the omega/analysis-moc-streamfunction branch from de20a34 to a48d53b Compare September 9, 2026 04:46
@xylar

xylar commented Sep 9, 2026

Copy link
Copy Markdown

All six of the previous review's findings are fixed — confirmed against the output rather than the diff, thank you. The time_bnds you were unsure about is right in a continuous run. It does have one problem, and chasing it turned up something more serious that is not about MOC at all: any stream carrying a temporal reduction segfaults when it writes into a file it has written before, which is what every restart causes. Details, a reproducer with MOC switched off entirely, and where I think the fault is, below.

Verified fixed

Re-ran the QU240 mock-up on this branch rebased onto omega/develop, built with intel (oneAPI 2025.2) on chrysalis, NumBins: 60, Regions: [Global], ReductionPeriod: [1Month], one simulated year. From ncdump -h on the first monthly file:

double MOCLatBinBoundaries(NMocLatBinBoundaries) ;
        MOCLatBinBoundaries:units = "degrees_north" ;
        MOCLatBinBoundaries:Units = "degrees_north" ;
double MOC_streamfunction_Global_TimeMean1Month(time, NumBinsLatCell_BinIndex, NVertLayersP1) ;
        MOC_streamfunction_Global_TimeMean1Month:units = "Sv" ;
        MOC_streamfunction_Global_TimeMean1Month:Units = "Sv" ;
        MOC_streamfunction_Global_TimeMean1Month:cell_methods = "time: mean" ;
        MOC_streamfunction_Global_TimeMean1Month:long_name = "Global Meridional Overturning Streamfunction" ;
double GeomZInterface_HorzMean(time, NVertLayersP1) ;
        GeomZInterface_HorzMean:units = "m" ;
        GeomZInterface_HorzMean:Units = "m" ;
double time(time) ;
        time:bounds = "time_bnds" ;
double time_bnds(time, D2) ;

The bin boundaries have lost their time dimension and now carry the operator's margin exactly — first boundary -90.000180000, last 90.000180000, matching (MaxLat - MinLat) * 1.0e-6 as CoordinateBinningOp applies it. Putting cell_methods in TimeMeanOp rather than in MOC.cpp was the right call: it covers every analysis group's time means.

In a continuous run time_bnds is exactly right. The first two monthly reductions of a run starting 0000-12-01 come out [0, 31] and [31, 62] days — correct spans, correctly abutting, sharing the time axis reference. Polaris reads all of this without special-casing anything, and the full analysis suite passes against a one-year mock-up.

One small note rather than a request: MOC.cpp repeats the margin formula with a comment pointing at CoordinateBinningOp.h instead of deriving it from the operator's own MinBin/BinWidth. The values agree today; the two expressions can still drift apart later.

Reduction streams segfault when they rewrite their own file

A continuation run crashes before its first time step. The trigger is not the restart itself but what the restart causes: the reduction alarm fires again at the restart instant and rewrites the output file the previous segment already wrote.

It is not specific to MOC. The tightest reproducer has the MOC group disabled entirely. Give GlobalStats a ReductionPeriod: [1Month], run a continuation from a restart twice in the same directory, and the second pass segfaults writing global_stats_1MonthTimeStats — the file the first pass created:

[info] [IOStream.cpp:2490] Successfully read stream RestartRead from file restart/ocn.restart.0001-02-01_00.00.00
[info] [IOStream.cpp:2806] Successfully wrote stream GlobalStats_1DayInstants to file global_stats_1DayInstants
srun: error: task 0: Segmentation fault      (exit 139)

The most informative line there is the one that succeeds. In the same pass, over the same restart, GlobalStats_1DayInstants rewrites its own pre-existing file without trouble. It is a snapshot stream, so it has no cell_methods and no time_bnds. The reduction stream beside it, writing into an equally pre-existing file, dies.

Seven runs, each continued from the identical 0001-02-01 restart:

build stream under test file already exists? result
this branch MOC, ReductionPeriod yes segfault
this branch MOC, ReductionPeriod no (fresh directory) completes
this branch MOC, SnapshotPeriod no completes
this branch MOC off, GlobalStats reduction, first pass no completes
this branch MOC off, GlobalStats reduction, second pass yes segfault
this branch, TimeDependent reverted to true MOC, ReductionPeriod yes segfault
previous tip of this branch, before the six commits MOC, ReductionPeriod yes completes

The last row is the attribution: the old code rewrites the same file in the same situation without crashing, so this is new in these commits. The sixth row rules out the TimeDependent=false change — I suspected that one first, and it is not the cause.

Where I think it is

The only new code in the write path is the time_bnds block in IOStream::writeStream, and it is also the only thing that distinguishes a reduction stream from a snapshot stream, since WriteTimeBnds is set precisely when some field carries cell_methods. That matches the evidence exactly: snapshot streams rewrite files fine, reduction streams do not.

The comment in that block states the assumption I think breaks:

// Define the CF-compliant time bounds variable (time_bnds) and attach the
// bounds attribute to the time variable. defineVar is called on every write
// (the file is re-entered in define mode each time) so TimeBndsID is valid
// for every frame; ...

When the file already exists, writeStream takes the readAllDims path rather than creating the file, so it is not in the same define state as a fresh write. If defineVar for an already-present time_bnds does not hand back a usable id under those conditions, TimeBndsID is stale or invalid by the time the data are written:

IO::writeNDVar(Bnds, OutFileID, TimeBndsID, Frame, &BndLengths);

which would put a bad variable id into the write and is consistent with a segfault rather than a clean error. I could not confirm this in a debugger — the compute nodes have no gdb, and eu-stack is blocked by ptrace restrictions — so treat the mechanism as indicated rather than proven. The evidence for where is strong; the exact failing call is inference.

Worth checking in the same pass: Dimension::create("D2", 2) is called before defineAllDims on every write, including writes to files that already have D2.

time_bnds after a restart is wrong even where it survives

The GlobalStats reduction that ran once and completed also shows what time_bnds does across a restart. It gives two records where one is expected:

record 0: time= 62.0 d   bnds=[  0.0,  62.0] d   span=62.0 d
record 1: time= 90.0 d   bnds=[ 62.0,  90.0] d   span=28.0 d

Record 1 is the February mean and is correct. Record 0 is written at the restart instant, duplicating an output time the previous segment already wrote, and claims to be a 62-day average of the whole run to date. Its value is not obviously wrong — KineticEnergyCell_SpatialMean is 0.00807 against the next month's 0.00842 — which makes it worse than an obvious one, because anything weighting reductions by the span in time_bnds will silently give it 62 days of weight.

The lower bound of 0 belongs to this PR. writeStream initialises PrevWriteTime from ModelClock->getStartTime() on a stream's first write, and neither PrevWriteTime nor FirstWrite survives a restart, so a continuation treats its first write as the first of the run. Clock::setCurrentTime(), which OceanInit calls when it reads a restart, resets the current time but deliberately leaves the clock's start time alone, so getStartTime() still returns the original simulation start.

Deriving the lower bound from the stream's own alarm interval rather than tracking it across writes — LowerBnd = ElapsedTimeR8 - <interval> — needs no state to be persisted, is restart-proof, and gives the same answer in a continuous run. That still leaves record 0 being written at all, which is the alarm firing at the restart instant and is not something this PR introduced.

A note on #524, which is not this PR's problem

While checking whether #524 might already address the restart behaviour, I tried merging it with this branch and hit conflicts in OceanInit.cpp, OceanDriver.h and omega_cxx2f_interface.cpp. They are not a disagreement between the two PRs. TimeInitParams and IOInitParams are already on develop, and this branch inherits them without touching either file; TimeStepperStartType exists only on omega/stop-time-changes, which is 179 commits behind develop. The conflict is between #524 and develop's current coupled-init signature. Flagging it only so it is not a surprise later.

How this was checked

The consumer is E3SM-Project/polaris#743, which reads this group's output, averages the reductions weighted by the length of each period, and plots the streamfunction. The full analysis suite passes against a one-year QU240 mock-up built from this branch. Work directories on chrysalis: the year-long run under /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260907/qu240-moc-rerun-20260909, the analysis under test_20260909/qu240-moc-analysis, and the restart experiments — each with a README.md saying what it isolates — under test_20260909/moc-restart-*.


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

@xylar

xylar commented Sep 13, 2026

Copy link
Copy Markdown

#548 adds a CFUnits class (src/infra/CFUnits.h) for deriving the units of a derived field: CFUnits::multiply("m s-1", "m2") gives m3 s-1, with the exponents of repeated symbols merged and unknown (empty) units propagating. The combineUnits in BinaryMultiplyOp here writes m s-1 * m2, which udunits does not parse, so cfchecks flags it; once #548 merges, the operators in this PR can use CFUnits instead. The same PR gives AnalysisOperator an inheritMetadata helper that copies units, standard_name and cell_methods from the input and appends the new reduction, which covers what BinnedAccumulatorOp, PrefixSumOp and ExtractRegionOp do by hand.


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

@brian-oneill

Copy link
Copy Markdown
Author

Fixed the restart crash by moving the Dimension::create earlier in writeStream before readAllDims populates AllDimIDs. Also made an edit that should fix the lower bound on time_bnds during the first write after a restart. There is still a spurious extra frame that gets written to the stream at the end of the first time step on a restart, which should be fixed by #524

@xylar xylar left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Ran a one-year QU240 run from d31f221 merged locally with #553 on Chrysalis (Intel, OpenMPI), at /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260914/qu240-monthly-avgs/ocean/analysis_test/QU240/forward. Both the MOC files and #553's monthly means now carry time_bnds with time:bounds pointing at it and the right interval, which is what #549 asks for, so this PR could say Fixes #549 in its description. One CF nit inline; the interval-end naming and time value are #554, not this PR's.


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

Comment thread components/omega/src/infra/IOStream.cpp Outdated
// Reuse the same units as the time variable (seconds since start).
std::string UnitString =
"seconds since " + StartTime.getString(4, 0, " ");
IO::writeMeta("units", UnitString, OutFileID, TimeBndsID);

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

cfchecks flags this: CF 7.1 says a boundary variable should not have a units attribute, since it takes the units of the variable it bounds. Dropping this line clears the warning.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

@brian-oneill, is this one you'd like to deal with here? I'll approve but I'd prefer having he units dropped here rather than having to do it in a follow-up PR. Bounds variables don't get units according to CF, they inherit them from the variable they are the bounds on.

@xylar

xylar commented Sep 15, 2026

Copy link
Copy Markdown

Testing: restart

I continued the one-year QU240 run (this branch at d31f221 merged with #553) from its 0001-03-01 restart for 61 days, on Chrysalis with Intel and OpenMPI, and compared the March and April files with the full run's:

/lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260915/qu240-restart-from-0001-03-01

With the clock started at the restart time, the monthly means and snapshots are bit for bit and time_bnds are the right intervals, [0, 31] and [31, 61] days. The April MOC is bit for bit; the March MOC, the first period after the restart, differs by 1.4e-14 in the streamfunction, so some state in the MOC chain is not the same after a restart as at a fresh start for its first period. Do you know what that would be?

With the clock left at the run's original start, every periodic alarm rang on the first step after the restart and the first time_bnds began at 0: Clock::setCurrentTime calls Alarm::updateStatus, which never moves RingTimePrev. #524 changes that call to Alarm::reset, so the time_bnds here are right after a restart only once #524 is in. Worth saying in the description. The files from that attempt are in attempt1_develop_clock/ there.


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

@brian-oneill

Copy link
Copy Markdown
Author

With the clock started at the restart time, the monthly means and snapshots are bit for bit and time_bnds are the right intervals, [0, 31] and [31, 61] days. The April MOC is bit for bit; the March MOC, the first period after the restart, differs by 1.4e-14 in the streamfunction, so some state in the MOC chain is not the same after a restart as at a fresh start for its first period. Do you know what that would be?

I replicated your test on Frontier with craygnu, and found everything BFB between analysis outputs of the continuous run and the restarts. The fact that the differences are near machine precision leads me to think it's Intel compiler optimizations that are responsible, though it's certainly strange it manifests specifically in the restart pipeline, and seems to wash out on subsequent outputs. Since the differences are so small and situationally specific, I lean toward tolerating the differences, while making a note of it.

@xylar

xylar commented Sep 25, 2026

Copy link
Copy Markdown

The fact that the differences are near machine precision leads me to think it's Intel compiler optimizations that are responsible, though it's certainly strange it manifests specifically in the restart pipeline, and seems to wash out on subsequent outputs. Since the differences are so small and situationally specific, I lean toward tolerating the differences, while making a note of it.

Okay, let's make a note of it. I think we won't typically rely on the analysis restart capability, so this may not matter anyway. We should also revisit after Phil's fix goes in because restarts are a little messy right now anyway.

@xylar xylar left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I'm approving based on my testing. This is excellent work!

See my final clean-up request above.

@xylar

xylar commented Sep 28, 2026

Copy link
Copy Markdown

I posted #578 to keep track of the restart issue.

@xylar

xylar commented Sep 28, 2026

Copy link
Copy Markdown

Testing

On Chrysalis with intel and OpenMPI, I tested the PR's test merge (c444c6e400) against a baseline built from develop at its merge base (c1aacdc2d7), using Polaris main (abe178ee8). All 50 CTests pass, and all 25 omega_pr tasks pass, including every baseline comparison.

Qualification: develop does not currently build standalone on Chrysalis (#572), so both the baseline and the PR build have #574 (edefec8790) merged in. The CTests also used the 260911 Icos480 sphere mesh from E3SM-Project/polaris#779, which develop has needed since #547.

Builds (-p):

baseline: /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260928/omega481/run-develop/build_omega/build_chrysalis_intel
PR:       /lcrc/group/e3sm/ac.xylar/polaris_1.1/chrysalis/test_20260928/omega481/run-pr481/build_omega/build_chrysalis_intel

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

@xylar

xylar commented Sep 28, 2026

Copy link
Copy Markdown

Aurora is down but tests on pm-cpu, pm-gpu, chrysalis, and frontier with CPU and GPU all passed (both CTests and omega_pr vs. develop as a baseline).

@xylar
xylar merged commit 36b1cd5 into E3SM-Project:develop Sep 28, 2026
8 checks passed
@xylar xylar mentioned this pull request Sep 29, 2026
7 tasks done
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.

3 participants