PerturbationAdvection: laterally-varying and FieldTimeSeries density#5596
Open
ewquon wants to merge 6 commits into
Open
PerturbationAdvection: laterally-varying and FieldTimeSeries density#5596ewquon wants to merge 6 commits into
PerturbationAdvection: laterally-varying and FieldTimeSeries density#5596ewquon wants to merge 6 commits into
Conversation
Per @jagoosw's review of #5587: - `to_intensive` / `to_extensive` now take `(ρ, ψ, i, j, k)` and own the field-indexing internally; each call site is a single expression rather than separately reading `ψ[iᴮ, jᴮ, kᴮ]` and passing the scalar through. - `step_right_open_boundary!` / `step_left_open_boundary!` drop the trailing `k` argument that duplicated `boundary_indices[3]`; the `step_right_boundary!` / `step_left_boundary!` aliases now forward through unchanged. - `_fill_*_halo!` callers updated to match the new signature. Pure refactor — behavior is unchanged on existing code paths (verified with the standard Bounded-grid PA smoke test). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Per @jagoosw's review of #5587: when the density field is at Centers and ψ is at a Face, the `ρψ → ψ → ρψ` conversion should use the density value at ψ's location, not the cell-center value half a step away. This commit adds a `_density_at_psi_location` helper that dispatches on ψ's instantiated location and applies the appropriate staggered-grid interpolator: - ψ at (Face, Center, Center) → ℑxᶠᵃᵃ - ψ at (Center, Face, Center) → ℑyᵃᶠᵃ - ψ at (Center, Center, Face) → ℑzᵃᵃᶠ - ψ at (Center, Center, Center) → direct ρ[i, j, k] `step_*_open_boundary!` and the `_fill_*_halo!` callers thread the ψ location through. The column-only `Field{Nothing, Nothing, Center}` case (anelastic reference profile) still works because the location-broadcasting makes the interpolated value identical to the column value. This loosens the previous hardcoded `ρ[1, 1, k]` indexing, which biased the boundary mass flux at every off-centerline face for laterally-varying density fields (e.g. fully-compressible model with a 3D prognostic ρ). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
The motivating use case for #5587 — regional atmospheric hindcasts driven by reanalysis — needs boundary density to vary in time and space, which a static `AbstractField` can't provide. This commit adds support for passing a `FieldTimeSeries` (or any `FlavorOfFTS`) as the `density` kwarg. The FTS conversion uses `interpolate(node(i, j, k, grid, loc...), Time(clock.time), ρ_fts, instantiated_location(ρ_fts), ρ_fts.grid)`, so the FTS may live on a different grid (and at a different location) than the simulation variable being converted. Time interpolation honors the FTS's `time_indexing` setting (`Clamp`/`Cyclical`/`Linear`). The FTS-specific overloads of `to_intensive` / `to_extensive` live in a new `src/OutputReaders/perturbation_advection_fts_density.jl` because `FlavorOfFTS`, `interpolate`, and `Time` aren't available at `BoundaryConditions` load time (OutputReaders is included afterwards). The static-`AbstractField` path defined in PA's file remains untouched. Docstring updated to document both density backends. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…variation
Covers three cases that the previous test suite didn't exercise:
- `test_perturbation_advection_constant_density`: uniformly-1 density given
as an `AbstractField`, a column-only `Field{Nothing, Nothing, Center}`
(anelastic reference profile), or a `FieldTimeSeries` must produce the
same boundary evolution as `density = nothing`. Validates the dispatch
paths for all three density types across the x, y, z orientations.
- `test_perturbation_advection_lateral_density`: a density Field that
varies along the boundary-parallel axis must produce a boundary
velocity that also varies along that axis. The pre-patch `ρ[1, 1, k]`
indexing would have produced a transverse-uniform boundary velocity;
this is a regression check against re-introducing that behavior.
Shared helpers `_pa_density_test_grid` and `_pa_density_step` factor out
the 3D grid + time-step setup. All wired into the existing
`@testset "Open boundary conditions"` alongside the FT loop.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
jagoosw
reviewed
May 18, 2026
| # Aliases for callers that follow the generic boundary-step naming. | ||
| # Default to a Center-located ψ so the existing column-density behavior | ||
| # (direct `ρ[i, j, k]` indexing) is preserved for any external caller. | ||
| @inline step_right_boundary!(bc::PAOBC, l, m, boundary_indices, boundary_adjacent_indices, |
Collaborator
There was a problem hiding this comment.
I'm not sure I understand the purpose of these methods?
Wouldn't the location never be (Center(), Center(), Center()) also for open boundaries?
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
2 participants
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.
Summary
Extends
PerturbationAdvection'sdensitykwarg to handle two cases the previous implementation could not:Laterally-varying density. The
ρψ → ψ → ρψconversion was hardcoded to readρ[1, 1, k]regardless of which boundary face was being processed, which biased the mass flux at every off-centerline face when ρ varied along the boundary. The conversion now reads ρ at ψ's grid location using standard staggered-grid interpolators (ℑxᶠᵃᵃ,ℑyᵃᶠᵃ,ℑzᵃᵃᶠ), so each face gets its own ρ value.FieldTimeSeriesdensity. Regional atmospheric hindcasts driven by reanalysis need boundary density to vary in time and space, which a staticAbstractFieldcouldn't supply. PA now accepts aFlavorOfFTSas thedensitykwarg, usinginterpolate(node(i, j, k, grid, loc...), Time(clock.time), ρ_fts, ...)so the density FTS may live on a different grid than the simulation and at a different location than ψ.Closes #5587.
Backwards compatibility
density = nothingpath is unchanged.Field{Nothing, Nothing, Center}case (anelastic reference profile, the documented use case from PR (0.105.1) Add density and gravity_wave_speed to PerturbationAdvection #5314) still works because location-broadcasting makes the staggered-grid interpolator return the same column value as directρ[1, 1, k]indexing did. Existing tests pass without modification.PerturbationAdvectionconstructor change.Implementation notes
The four commits are reviewable independently:
PerturbationAdvection: support 3D / laterally-varying density at open boundaries #5587):to_intensive/to_extensivenow own the ψ-indexing;step_*_open_boundary!drop a spurious trailingkargument that duplicatedboundary_indices[3]. Pure refactor — behavior unchanged._pa_density_valuedispatched oninstantiated_location(ψ). Threads the ψ location through the kernel chain.FlavorOfFTSoverload ofto_intensive/to_extensivein a newsrc/OutputReaders/perturbation_advection_fts_density.jl(kept there becauseFlavorOfFTS,interpolate, andTimearen't available atBoundaryConditionsload time).Test plan
Pkg.test("Oceananigans"; test_args=["--group=time_stepping_2"])passes locally on CPU (333/333 in the boundary-conditions integration testset).test_perturbation_advection_constant_densityvalidates that constant ρ=1 given as anAbstractField, a column-onlyField{Nothing, Nothing, Center}, or aFieldTimeSeriesproduces the same boundary evolution asdensity = nothing, across all three orientations (x, y, z).test_perturbation_advection_lateral_densityis a regression check: a density field varying along the boundary-parallel axis must produce a boundary velocity that also varies — pre-patch code would have produced a transverse-uniform velocity.density = nothing).