Add Omega support for the seamount tests - #678
Open
xylar wants to merge 14 commits into
Open
Conversation
Collaborator
Author
|
@cbegeman, feel free to have a quick look but I need to look this over carefully myself and don't have time today. I just wanted to get a draft PR up to indicate that the work was mostly done so you didn't also try to do it. |
xylar
force-pushed
the
add-seamount-omega-support
branch
from
July 30, 2026 16:54
71828c0 to
b3bbdde
Compare
xylar
commented
Jul 30, 2026
cbegeman
reviewed
Jul 30, 2026
cbegeman
reviewed
Jul 30, 2026
5 tasks
The seamount initial condition back-solved temperature from the target
Beckmann and Haidvogel density using only
rho = rho_ref - alpha * (T - T_ref)
dropping the beta * (S - S_ref) term that compute_linear_density() and
both ocean models actually apply. With beta = 0.8 and S = 35 PSU the
density the model saw was 28 kg m-3 above the profile the test intended
to impose. Under MPAS-Ocean that is a constant offset and harmless, but
Omega's p-star coordinate maps geometric to pseudo thickness using the
absolute density, so it is not harmless there.
Separately, seamount.cfg set eos_linear_Tref = 5, which Omega rejects:
its linear EOS has no reference temperature or salinity, and
OceanModelStep._get_linear_eos_replacements() raises for a nonzero Tref
or Sref. This blocked Omega support outright.
Fold the reference state into rhoref at Tref = Sref = 0
(1028.0 + 0.2 * 5.0 - 0.8 * 35.0 = 1001.0) and back-solve through the
full linear EOS. The two changes cancel exactly in the temperature
field, so temperature, salinity and velocity are unchanged to machine
precision; only the model's density diagnostic moves, by -28 kg m-3,
onto the intended Beckmann and Haidvogel profile.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Follow the overflow pattern of a subdirectory per vertical coordinate, giving planar/seamount/sigma and planar/seamount/zstar. Sigma remains the main target; z-star is the control, since the flatter its layers the smaller the spurious pressure-gradient velocity should be. Both are geometric coordinates, so both will serve Omega once the forward step supports it. coord_type moves out of the shared seamount.cfg into per-coordinate seamount_sigma.cfg and seamount_zstar.cfg. The z-star tree uses partial bottom cells. Sigma ignores partial cells, but z-star cuts the reference grid at the seafloor, and on this bathymetry that needs snapping to be usable: sampling the seamount profile densely at 32 levels gives a minimum bottom-cell thickness of 7 cm with partial_cell_type = None, against 15.6 m with "partial" and min_pc_fraction = 0.1. min_pc_fraction is added to the shared config, following the convention of the other ocean tasks. Also adopt the usual os.path.join(component.name, ...) form for the shared config filepath, which the seamount had been omitting. No test suite references the seamount, so the task path change is safe. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Split forward.yaml into shared "ocean:", "mpas-ocean:" and "Omega:" sections following the overflow task, and teach forward.py to configure either model. Omega has no split-explicit integrator, so it is limited by the barotropic gravity wave: sqrt(g * 5000 m) = 221 m/s against dx = 6.7 km gives a CFL-1 limit near 4.5 s/km, where MPAS-Ocean's split-explicit step runs at 10. The models therefore need separate time steps and integrators, so dt_per_km/time_integrator gain omega_ counterparts. The Omega step starts at 2.5 s/km, matching MPAS-Ocean's own barotropic step; this is a first guess pending a stability test. Vertical mixing is switched off. The exact solution is a resting ocean, so the only thing that would ever trigger convection here is a spurious pressure gradient -- the quantity being measured -- and convection is a known source of MPAS-Ocean/Omega divergence, since Omega uses a single coefficient for both convective diffusivity and viscosity. cvmix stays enabled with its coefficients zeroed rather than disabled, because config_use_cvmix = false falls back on MPAS-Ocean's constant vertical viscosity and diffusivity, which are not zero. This changes MPAS-Ocean answers. Bottom drag stays implicit and constant in both models, which the recently added Omega implicit bottom drag support makes possible, and moves to the shared section. Horizontal tracer advection order is pinned so the two models do not use their differing defaults. Verified with polaris setup for Omega: the generated omega.yml has TimeStepper RungeKutta4, TimeStep 16.750 s, BottomDragTendency implicit and constant, VertMix background zeroed with convection and shear off, and Eos linear with RhoT0S0 1001.0. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Both changes re-baseline MPAS-Ocean. Cd 1e-2 -> 1e-3. The implicit drag damping timescale is h_bot / (Cd * |u|), and with sigma at 32 levels the bottom layer over the seamount summit -- where the coordinate is most tilted and the signal lives -- is only ~16 m thick. At 1e-2 a 5 cm/s spurious velocity damps in 0.4 days against a 6 day run, while a 1 cm/s one damps in 1.8 days; the metric saturates, capping a bad pressure gradient while leaving a good one untouched. At 1e-3 drag remains a backstop against unbounded spin-up without flattening the signal. Comparable idealized planar cases in Compass use 1e-2, but none of them is measuring a quantity that drag suppresses. vert_levels 10 -> 32. A multiple of 16, which Omega prefers, dividing the 5000 m bottom depth into exact 156.25 m layers, and in the range the literature uses for this case (20 in Beckmann and Haidvogel, 20-30 in Shchepetkin and McWilliams) rather than below it. It also rescues the z-star variant: at 10 uniform levels the 500 m summit column has a single level, against three full levels plus a partial at 32. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The seamount hands Omega a geometric coordinate and relies on Polaris converting it to pseudo-height, rather than on any p-star-specific initialization. These tests exercise that conversion directly, on a seamount-like set of columns, with no file I/O or step framework. They cover: - the temperature back-solve reproduces the Beckmann and Haidvogel density through compute_linear_density(); - geometric -> pseudo -> geometric recovers the column thickness, for sigma and for z-star with partial cells; - sigma layers are exactly proportional to column depth; - interior interfaces sit at genuinely different pressures from column to column, which is the property a clipped reference grid would destroy by leaving the whole pressure-gradient signal in the bottom cell. The last of these gets an explicit negative control on a flat bottom, where the same measurement must come back at round-off. Without it the tilt assertion could not fail, and a guard that cannot fire is worse than no guard. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Replace "These tasks support only MPAS-Ocean", describe the two coordinate trees, and explain the geometric-to-pseudo-height route that lets one init step serve both models. Record the decisions that are not obvious from the code: why eos_linear_Tref must be zero and the reference state is folded into eos_linear_rhoref; why the temperature back-solve has to include the salinity term; why all vertical mixing is off but cvmix stays enabled; why the bottom drag coefficient sits an order of magnitude below the value comparable idealized planar cases use; why z-star needs partial cells; and why Omega takes a shorter time step. Also drop config options from the User's Guide that the task does not have: seamount_density_ref, seamount_density_Tref and seamount_density_alpha, replaced by the eos_linear_* options in the ocean section, and coriolis_parameter, now type and constant_f in the coriolis section. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add `ct_from_potential_density()` to `polaris/ocean/eos/teos10.py`, the TEOS-10 counterpart of back-solving a linear equation of state for temperature. It lets a test case specify its stratification as a density profile rather than a temperature profile, which several idealized cases do -- the seamount's Beckmann and Haidvogel profile among them. Referencing the target density to the surface makes the inversion exact and non-iterative: `gsw.CT_from_rho` is called at `p = 0`, so the result does not have to be folded into the pressure/specific-volume iteration that the geometric-to-pseudo-height conversion already performs. Only the cold branch of `CT_from_rho` is returned. The warm branch exists only below the density maximum, which for seawater lies below the freezing point at oceanographic salinities; `gsw` reports it as NaN throughout the range any of these cases use. An unattainable target density raises rather than returning NaN, since a NaN would otherwise propagate silently into an initial condition. NaN input, which marks below-bottom cells, still passes through untouched. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Split the seamount into `planar/seamount/{linear,nonlinear}/{sigma,zstar}`,
following the `planar/overflow/{eos}/{coord}` pattern. The nonlinear trees
use TEOS-10 for Omega and Jackett-McDougall for MPAS-Ocean, selected by the
shared `polaris.ocean.eos` `teos10.cfg`; no new EOS or vertical-coordinate
plumbing is needed, because the geometric-to-pseudo-height conversion that
Omega's coordinate goes through already dispatches on `eos_type`.
The Beckmann and Haidvogel initial condition is specified as a density
profile, and the linear EOS is inverted for the temperature that reproduces
it. That inverse does not carry over: TEOS-10 in-situ density at 5000 m is
near 1050 kg m^-3 from compression alone, so the 1025-1028 kg m^-3 profile
cannot be read as in-situ density at any temperature. It is therefore read
as a potential density referenced to the surface and inverted with
`ct_from_potential_density()`. Referencing to the surface is not just a
convenience: it makes the linear and nonlinear trees share a buoyancy
stratification exactly, so a difference in spurious velocity between them is
attributable to the equation of state rather than to a different N^2.
Verified on the 32-level grid: sigma_0 agrees between the two trees to
2.3e-13 kg m^-3 while in-situ density spans 1025.1-1050.6, and the geometric
layer positions are identical to the last bit.
The profile and tracer construction move to a new `init_utils.py`, a leaf
module free of `mpas_tools` and of the step framework so the unit tests can
import it, named after overflow's module of the same name. The linear EOS
coefficients move from the shared `seamount.cfg` into a new
`seamount_linear.cfg`, since they mean nothing to the nonlinear trees.
For MPAS-Ocean the CT/SA tracers have to become potential temperature and
practical salinity. The step does not do that itself: it leaves the tracers
in the convention `eos_type` implies and `write_initial_state_dataset()`
converts them at the `[ocean]` nominal lon/lat, computing the pressure the
salinity correction needs from the layer thicknesses, since a geometric
vertical coordinate carries none.
The tracers now carry the `Time` dimension that `zMid`, `layerThickness` and
the rest of the state already had, matching overflow. TEOS-10 requires its
inputs to be aligned and the linear EOS did not, which is why this went
unnoticed. Values are unchanged: `init.nc` and `vert_coord.nc` for both
linear trees match the test_20260730 baseline exactly across all 23 and 5
variables respectively, with only the added dimension differing.
Paths change for the existing trees, from `planar/seamount/{sigma,zstar}` to
`planar/seamount/linear/{sigma,zstar}`. No suite references the seamount.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Cover both equation-of-state branches of the seamount tracer construction, in the self-contained style of the existing coordinate tests: a synthetic two-column set spanning the seamount's depth range, no file I/O and no step framework. The tests that matter are the two that state the design: - both branches reproduce the same Beckmann and Haidvogel profile, the linear one as density and the nonlinear one as potential density at zero reference pressure, agreeing to 1e-11 kg m^-3. This is the property that makes a difference in spurious velocity between the trees attributable to the equation of state rather than to a different N^2. - the negative control: in-situ density under TEOS-10 departs from the surface-referenced profile by more than 20 kg m^-3 at 5000 m, monotonically with depth, against the 3 kg m^-3 the whole profile spans. Without that departure the nonlinear tree would measure nothing the linear one does not. Also covers the linear back-solve against `compute_linear_density()`, which the coordinate tests did not reach, the linear stratification branch, and both unsupported-option errors. The conversion of the nonlinear tracers to the MPAS-Ocean convention is the framework's and is covered by the framework's own tests. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Replace the two-tree coordinate description with the four-tree
`{linear,nonlinear}/{sigma,zstar}` layout, and split the "variants" section
into what the coordinate controls (how much the layers tilt) and what the
equation of state controls (whether density depends on pressure).
The initial-conditions section now explains why the Beckmann and Haidvogel
profile is read as a potential density referenced to the surface under a
nonlinear equation of state -- in-situ density at 5000 m is near
1050 kg m^-3, so no temperature reproduces the 1025-1028 profile below about
1000 m -- and what the surface reference buys: an identical buoyancy
stratification between the linear and nonlinear trees.
Also records the caveat that TEOS-10 and Jackett-McDougall are different
functions rather than two implementations of one, so the two models are
only expected to agree qualitatively in the nonlinear trees.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
xylar
force-pushed
the
add-seamount-omega-support
branch
from
August 1, 2026 11:22
b3bbdde to
5ffd34e
Compare
Omega has no split time stepper, so MPAS-Ocean gives up its split-explicit one and both models now run RK4 at the same `dt_per_km = 4.0`. The per-model `omega_dt_per_km` and `omega_time_integrator` config options are removed; only the integrator name is still translated for Omega (`RK4` -> `RungeKutta4`). `btr_dt_per_km` is kept but is inert unless `time_integrator` is set back to a split-explicit scheme. Vertical mixing is also switched off for MPAS-Ocean with `config_use_cvmix = false`, matching Omega, whose coefficients were already zeroed. The comment claiming that disabling cvmix falls back on a nonzero constant vertical viscosity and diffusivity is removed: in the current MPAS-Ocean, `ocn_vmix_coefs` zeroes the coefficients and `ocn_vmix_coefs_cvmix_build` returns immediately when cvmix is off, so nothing writes into them. The implicit vertical mixing solve itself has to stay enabled in both models, since that is what applies the implicit bottom drag, which is now noted where it would otherwise be tempting to disable it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The default seamount task integrates for 6 days, which is long enough that it is not worth making the user re-run the task just to get the plots. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The default seamount task runs for 6 days, which is far too long for a pull request suite. Add a `short` task alongside it in each of the four trees: the same mesh, initial condition, time step and physics, run for 1 hour. An hour is too short for the spurious circulation to develop, so it measures nothing about the pressure gradient error, but it is enough to catch a change in the answer cheaply. `viz` does not run by default in the short task, unlike the default task, since re-running an hour-long forward step to get the plots costs almost nothing. Two of the four are added to `mpaso_pr` and `omega_pr`, `linear/zstar` and `nonlinear/sigma`, which covers both equations of state and both vertical coordinates in two runs rather than four. `run_duration` is now given in hours rather than days so that both task sections can use the same units as `output_interval`. The default task is unchanged at 144 hours. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
6.7 km was never prescribed by Beckmann and Haidvogel; it was carried over from the legacy Compass version of this test, where it appears to have been chosen as roughly 320 / 48. It does not divide the domain, and since the mesh is built from `nx * resolution` rather than from `lx`, the domain was really 321.6 km by 324.9 km rather than the 320 km by 320 km the config and docs advertise. 320 km divides into exactly 50 cells at 6.4 km, so `lx` now means what it says (`ly` is still approximate at 321.5 km, since the hexagonal mesh spacing in y is sqrt(3)/2 times the spacing in x). The test itself is unchanged in any way that matters. What the pressure gradient error responds to is the discrete steepness of the seamount, max |dh| / (h1 + h2) between adjacent cells, which for a Gaussian seamount is about 1.5 * resolution / seamount_width: 0.22 at 6.7 km and 0.21 at 6.4 km. Both bracket the value quoted for the configuration in the literature, so the case is neither easier nor harder in any meaningful sense, and the bottom drag and time step reasoning carries over untouched. `dt_per_km` in particular needs no retuning: the gravity wave frequency goes like 1 / dcEdge, so the stability cap of 5.2 s/km is the same at any resolution. The cost is 8% more cells (2900 rather than 2688) and a 4% shorter time step. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Collaborator
Author
xylar
marked this pull request as ready for review
August 1, 2026 15:07
Collaborator
Author
|
@cbegeman, this is ready to review. It's not urgent but seems nice to have! |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.








The
seamounttask measures the spurious velocity produced by the pressure gradient error in tilted layers. This PR makes it run under Omega as well as MPAS-Ocean, splits it into a tree per equation of state and per vertical coordinate followingoverflow, and revisits the physics and numerics it inherited from Compass so that the signal it measures is not suppressed and the two models are comparable.Task trees
planar/seamount/{sigma,zstar}toplanar/seamount/{linear,nonlinear}/{sigma,zstar}, following theplanar/overflow/{eos}/{coord}pattern. No suite referenced the seamount before this PR.sigmais the main target (every layer tilted by construction);zstaris the control, where layers are nearly level and the spurious velocity should be much smaller.linearhas no pressure dependence;nonlinearuses TEOS-10 for Omega and Jackett-McDougall for MPAS-Ocean. The two are different functions, not two implementations of one, so the models are only expected to agree qualitatively in the nonlinear trees.coord_type/partial_cell_typemove toseamount_{sigma,zstar}.cfgand theeos_linear_*options toseamount_linear.cfg, layered over the sharedseamount.cfgalong withlinear.cfgorteos10.cfgfrompolaris.ocean.eos.partial_cell_type = Noneagainst 15.6 m withpartialandmin_pc_fraction = 0.1.Omega support
initstep serves both models: the framework convertsrestingThickness/layerThicknesstoRefPseudoThickness/PseudoThicknessby integrating gauge pressure down each column. No p-star-specific init and no reference pseudo-grid, soPStarInitStepis not used.forward.yamlsplits into sharedocean:,mpas-ocean:andOmega:sections;forward.pyconfigures either model.dt_per_km = 4.0, against a barotropic-gravity-wave stability cap of 5.2 s/km.btr_dt_per_kmis kept but inert unlesstime_integratoris set back to a split-explicit scheme.SpecVolis required in the OmegaHistorystream, since Polaris reads back a geometriclayerThicknessasRhoSw * SpecVol * PseudoThickness.Initial condition
beta * (S - S_ref)term thatcompute_linear_density()and both models actually apply, leaving the model's density 28 kg m^-3 above the intended Beckmann and Haidvogel profile — harmless in MPAS-Ocean, not in Omega, where the pseudo-height mapping depends on absolute density. It now inverts the full linear EOS.eos_linear_Tref = 5is rejected outright for Omega, whose linear EOS has no reference temperature or salinity, so the reference state is folded intorhorefatTref = Sref = 0(1028.0 + 0.2 * 5.0 - 0.8 * 35.0 = 1001.0). This cancels exactly against the fix above: temperature, salinity and velocity are unchanged to machine precision and only the density diagnostic moves, onto the intended profile.ct_from_potential_density()inpolaris/ocean/eos/teos10.py— exact and non-iterative atp = 0, cold branch only, raising rather than returning NaN for an unattainable density.sigma_0agrees to 2.3e-13 kg m^-3 on the 32-level grid) while in-situ densities differ by more than 20 kg m^-3 at depth, so a difference in spurious velocity between the trees is attributable to the EOS rather than a differentN^2.init_utils.py, a leaf module free ofmpas_toolsand the step framework so the unit tests can import it. The step leaves tracers in the conventioneos_typeimplies;write_initial_state_dataset()converts CT/SA to PT/practical salinity for MPAS-Ocean, computing pressure from layer thicknesses since a geometric coordinate carries none.Timedimension the rest of the state already had, which TEOS-10 requires for alignment and the linear EOS did not.update_namelist_eos()already mapsteos-10tojm, and the pseudo-height conversion already dispatches oneos_type.Physics and resolution
1e-2→1e-3: the damping timescale ish_bot / (Cd * |u|)and the sigma bottom layer over the summit is only ~16 m thick, so at1e-2a 5 cm/s spurious velocity damps in 0.4 days against a 6 day run — capping a bad pressure gradient while leaving a good one untouched. Comparable Compass cases use1e-2, but none measures a quantity drag suppresses.nx * resolution, the domain was really 321.6 km across, not the advertised 320 km. Discrete seamount steepness, which is what the pressure gradient error responds to, barely moves (0.22 → 0.21). Costs 8% more cells and a 4% shorter time step.config_use_cvmix = falseto match Omega. The implicit vertical mixing solve itself stays on in both, since that is what applies the implicit bottom drag.Short task and suites
shorttask is added in each tree: same mesh, initial condition, time step and physics, run for 1 hour. Too short to develop the spurious circulation, but enough to catch a change in the answer cheaply.linear/zstar/shortandnonlinear/sigma/shortare added tompaso_prandomega_pr, covering both equations of state and both coordinates in two runs. Nothing is in the nightly suites yet, since thedefaulttasks are 6 day runs.run_durationis now in hours rather than days, matchingoutput_interval; the default task is unchanged at 144 hours.viznow runs by default indefault(a 6 day run is not worth repeating for plots) and stays opt-in inshort.Tests
Three new self-contained unit test modules (no file I/O, no step framework), 19 tests:
tests/ocean/eos/test_ct_from_potential_density.py— round trips throughgsw.rhoatp = 0, the unattainable-density error, NaN pass-through.tests/ocean/seamount/test_sigma_to_pseudo.py— geometric round trip for sigma and z-star with partial cells, sigma layers proportional to column depth, and interior interfaces at genuinely different pressures from column to column (the property a clipped reference grid would destroy), with a flat-bottom negative control so the tilt assertion can fail.tests/ocean/seamount/test_tracers.py— both EOS branches reproduce the same profile to 1e-11 kg m^-3; negative control that TEOS-10 in-situ density departs from it by more than 20 kg m^-3 at 5000 m against the 3 kg m^-3 the profile spans; the linear back-solve againstcompute_linear_density(); both unsupported-option errors.Answer changes
MPAS-Ocean answers change from the vertical refinement, the resolution change, the lower bottom drag, RK4 in place of split-explicit, and disabling cvmix.
init.ncandvert_coord.ncfor both linear trees match thetest_20260730baseline exactly (23 and 5 variables), with only the addedTimedimension differing.Checklist
api.md) has any new or modified class, method and/or functions listedTestingcomment in the PR documents testing used to verify the changes