Skip to content

Terrain-following: sample topography at ξnode/ηnode for lat-lon support#776

Closed
ewquon wants to merge 3 commits into
mainfrom
eq/fix-terrain-following-latlon-coords
Closed

Terrain-following: sample topography at ξnode/ηnode for lat-lon support#776
ewquon wants to merge 3 commits into
mainfrom
eq/fix-terrain-following-latlon-coords

Conversation

@ewquon

@ewquon ewquon commented Jun 10, 2026

Copy link
Copy Markdown
Collaborator

Summary

set_topography! sampled the user terrain function at xnode(i, grid, Center()) / ynode(j, grid, Center()), which are RectilinearGrid-oriented and don't generalize: on a LatitudeLongitudeGrid those accessors don't hand the topography function the horizontal (λ, φ) it expects.

This switches to ξnode/ηnode — the grid-agnostic horizontal node accessors:

  • (x, y) on a RectilinearGrid
  • (λ, φ) on a LatitudeLongitudeGrid → terrain-following now works on lat-lon LAMs
  • nothing for the degenerate meridional node on a 2D x–z (Periodic, Flat, Bounded) grid — exactly matching prior behavior

Why not node(...)[1:2]?

The obvious grid-agnostic substitute is wrong on Flat grids: node drops Flat directions from its returned tuple, so on a 2D x–z grid it yields (x, z) and leaks the vertical coordinate into the meridional slot. ξnode/ηnode are the primitives node is built from but never drop Flat directions. A regression test pins this down.

Tests

  • Horizontal-coordinate handling across Flat / 3D / lat-lon grids. The Flat-grid case uses a probe h(x, y) = y === nothing ? 1.0 : 2.0 asserting the second argument is nothing (not z) — this assertion fails against the node[1:2] approach.
  • LatitudeLongitudeGrid terrain-physics testset mirroring the existing rectilinear one: follow_terrain! + CompressibleDynamics(...; terrain_metrics) + one time_step! → finite w and contravariant vertical velocity. Cheap insurance that the metric terms behave under spherical geometry.

Full test/terrain_following.jl passes 664/664; quality_assurance (ExplicitImports/Aqua) passes 36/36.

🤖 Generated with Claude Code

`set_topography!` sampled the user topography at `xnode(i, grid, Center())`
/ `ynode(j, grid, Center())`, which are RectilinearGrid-oriented and don't
generalize: on a LatitudeLongitudeGrid those accessors don't give the
horizontal (λ, φ) the topography function expects.

Sample at `ξnode`/`ηnode` instead — the grid-agnostic horizontal node
accessors. They return `(x, y)` on a RectilinearGrid and `(λ, φ)` on a
LatitudeLongitudeGrid, so terrain-following now works on lat-lon LAMs. On
a 2D x–z (Periodic, Flat, Bounded) grid the second argument stays the
degenerate meridional node `nothing`, exactly matching the previous
behavior.

Note `node(...)[1:2]` is *not* a correct grid-agnostic substitute: `node`
drops Flat directions from its tuple, so on a 2D x–z grid it would yield
`(x, z)` and leak the vertical coordinate into the meridional slot.
`ξnode`/`ηnode` are the primitives `node` is built from but never drop
Flat directions; a regression test pins this down.

Tests:
- horizontal-coordinate handling across Flat / 3D / lat-lon grids, with a
  Flat-grid probe asserting the second argument is `nothing` (not z)
- a LatitudeLongitudeGrid analogue of the terrain-physics testset
  (follow_terrain! + CompressibleDynamics + one time_step! → finite w),
  exercising the metric terms under spherical geometry

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@@ -98,7 +98,12 @@ end
# is intentionally more general.
function set_topography!(h_field, grid, topography::Function)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

what is wrong with using set!(h_field, topography) here? Is there a difference between h_field.grid and the argument grid?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

This is more robust, per the comment below on Flat dims. Good call on the redundant grid arg

@ewquon ewquon Jun 10, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Ah, I missed the ::Function qualifier. I've added a generic version, which is actually what would be called when using real terrain data:

regrid_topography(grid; ETOPO2022()) → field → follow_terrain!(grid, field) → set_topography! → set!

@codecov

codecov Bot commented Jun 10, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

ewquon and others added 2 commits June 9, 2026 23:29
Per review feedback on set_topography!:

- `h_field.grid === grid` (the field is built from `grid`), so the `grid`
  parameter was redundant — derive it from `h_field.grid` and drop the arg.

- The old comment claimed `set!` "evaluates on-device, which would fail for
  non-GPU-compatible user functions". That is stale: Oceananigans'
  `set_to_function!` already evaluates on a CPU grid and transfers to the
  device. The real reason to keep the manual sampling is that `set!` evaluates
  a FunctionField by splatting `node(i,j,k,grid,…)` into the user function,
  which ties arity to the grid dimensionality (`f(x, y, z)` on a 3D grid) and,
  on a 2D x–z Flat grid, drops the Flat direction and calls `f(x, z)` — leaking
  z into the meridional slot. Sampling ξnode/ηnode keeps the `h(x, y)` API and
  always passes horizontal-only coordinates. Comment rewritten accordingly.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The `::Function` method needs manual ξnode/ηnode sampling to keep the
`h(x, y)` API and avoid `set!`'s grid-dimension-dependent arity. But array,
Field, and number topographies need none of that — `set!` copies/broadcasts
them onto `h_field` directly. Add a generic fallback that delegates to `set!`,
so `follow_terrain!` accepts precomputed height fields and constants, not just
functions.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@ewquon ewquon requested a review from glwagner June 10, 2026 05:50
@glwagner

Copy link
Copy Markdown
Member

This was superceded by #712 (the work here ended up being incorporated into #712, which was easier than merging this PR separately). Thank you @ewquon and @kaiyuan-cheng for your work on this!

@glwagner glwagner closed this Jun 14, 2026
@glwagner glwagner deleted the eq/fix-terrain-following-latlon-coords branch June 14, 2026 06:04
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.

2 participants