Skip to content

Add FieldTimeSeriesRelaxation#562

Open
ewquon wants to merge 27 commits into
mainfrom
eq/newton_relax
Open

Add FieldTimeSeriesRelaxation#562
ewquon wants to merge 27 commits into
mainfrom
eq/newton_relax

Conversation

@ewquon
Copy link
Copy Markdown
Collaborator

@ewquon ewquon commented Mar 14, 2026

This implements a Newtonian relaxation forcing term, to be used for grid/analysis nudging. The example script depends on ERA5 data retrieval developed in NumericalEarth/NumericalEarth.jl#93.

In the single-column example, nudging is applied above 1.5km to a domain with slip-wall boundary and uniform initial fields.

6-hour timescale:
single_column_nudging_tau6h

1-hour timescale (more responsive):
single_column_nudging_tau1h

ewquon and others added 6 commits March 13, 2026 19:06
Implements a generalized nudging forcing that relaxes any prognostic
field toward a reference FieldTimeSeries with a configurable time scale
and cutoff height below which no nudging is applied.

Two modes are supported: profile mode (compares the domain horizontal
average to a reference column, for LES grid nudging) and 3D mode
(compares each local value to the reference, for mesoscale nudging).
The mode is selected automatically based on whether the reference FTS
has horizontal dimensions and whether a reference_position is provided.

Also adds `clock` to the forcing materialization context, which is
needed by RelaxationForcing (and available for future forcing types)
to access the current simulation time in compute_forcing!.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Tests cover stub construction, profile and 3D mode materialization,
analytical verification of the nudging tendency magnitude, and the
z_bottom cutoff (zero tendency below the cutoff height).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
The clock was added to forcing_context but not to the function
signatures or call site, causing an UndefVarError at model construction.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Demonstrates RelaxationForcing for u, v, θ above 1500 m toward hourly
ERA5-style reference FieldTimeSeries, combined with time-varying surface
heat and moisture flux boundary conditions from scalar FieldTimeSeries.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Comment thread src/Forcings/relaxation_forcing.jl Outdated
Comment thread src/Forcings/relaxation_forcing.jl Outdated
Comment thread src/Forcings/relaxation_forcing.jl Outdated
Comment thread examples/single_column_nudging.jl Outdated
Comment thread examples/single_column_nudging.jl Outdated
Comment thread examples/single_column_nudging.jl Outdated
Comment on lines +79 to +83
Δλ = temp.grid.Δλᶜᵃᵃ
Δφ = temp.grid.Δφᵃᶜᵃ

xn = [ref_loc.longitude - Δλ/2, ref_loc.longitude + Δλ/2]
yn = [ref_loc.latitude - Δφ/2, ref_loc.latitude + Δφ/2]
Copy link
Copy Markdown
Member

@navidcy navidcy Mar 14, 2026

Choose a reason for hiding this comment

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

we can use xnodes and ynodes as well?
or just

xn, yn, zn = Oceananigans.nodes(temp.grid, Center, Center, Face, with_halos=false)

?

Comment thread examples/single_column_nudging.jl
Comment thread examples/single_column_nudging.jl Outdated
@codecov
Copy link
Copy Markdown

codecov Bot commented Mar 14, 2026

Codecov Report

❌ Patch coverage is 58.82353% with 35 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
src/Forcings/fts_relaxation.jl 58.44% 32 Missing ⚠️
...c/BoundaryConditions/thermodynamic_variable_bcs.jl 50.00% 3 Missing ⚠️

📢 Thoughts on this report? Let us know!

Comment thread src/Forcings/relaxation_forcing.jl Outdated
@glwagner
Copy link
Copy Markdown
Member

I suggest we call this HorizontalAverageRelaxation. In a future PR we can consider merging with Oceananigans.Forcing.Relaxation.

@ewquon
Copy link
Copy Markdown
Collaborator Author

ewquon commented Mar 17, 2026

I suggest we call this HorizontalAverageRelaxation. In a future PR we can consider merging with Oceananigans.Forcing.Relaxation.

@glwagner This naming only captures half the functionality, for an idealized simulation configuration. The more general use case of this forcing is to nudge a larger (mesoscale) domain with open boundaries towards a (possibly coarser) background field to ensure that the solution remains physically consistent with the larger-scale flow as it evolves in time and space.

How about GridNudging which is consistent with meteorological community usage?

@glwagner
Copy link
Copy Markdown
Member

I suggest we call this HorizontalAverageRelaxation. In a future PR we can consider merging with Oceananigans.Forcing.Relaxation.

@glwagner This naming only captures half the functionality, for an idealized simulation configuration. The more general use case of this forcing is to nudge a larger (mesoscale) domain with open boundaries towards a (possibly coarser) background field to ensure that the solution remains physically consistent with the larger-scale flow as it evolves in time and space.

How about GridNudging which is consistent with meteorological community usage?

That works for me. However, I do not think we should implement the more general 3D nudging in this feature. Instead, I think we should merge this feature with Oceananigans.Forcing.Relaxation first. After we merge the two, we can implement the 3D nudging to a 3D FieldTimeSeries in Oceananigans. It will also be useful for open boundary conditions for ocean models.

What do you think @ewquon ? The point of proposing HorizontalAverageRelaxation was simply to propose a temporary name which would be soon deleted when this is merged into Oceananigans.

@ewquon
Copy link
Copy Markdown
Collaborator Author

ewquon commented Mar 17, 2026

Well the thing is @glwagner, the 3D nudging has already been implemented here (it just needs more testing)...

# 3D mode (reference_column = Nothing): compare local value to 3D reference
@inline function (f::RelaxationForcing{<:Any, <:Any, <:Any, <:Any, <:Any, Nothing})(i, j, k, grid, clock, fields)
z = znode(i, j, k, grid, Center(), Center(), Center())
ϕ = @inbounds f.current_field[i, j, k]
ϕᵣ = @inbounds f.target[i, j, k]
ρ = @inbounds f.density[1, 1, k]
tendency = -ρ *- ϕᵣ) / f.time_scale
return ifelse(z < f.z_bottom, 0, tendency)
end
# Profile mode (reference_column = NTuple{2,Int}): compare horizontal average to reference column
@inline function (f::RelaxationForcing{<:Any, <:Any, <:Any, <:Any, <:Any, <:NTuple{2}})(i, j, k, grid, clock, fields)
z = znode(i, j, k, grid, Center(), Center(), Center())
ϕ̄ = @inbounds f.current_field[1, 1, k]
ϕᵣ = @inbounds f.target[1, 1, k]
ρ = @inbounds f.density[1, 1, k]
tendency = -ρ * (ϕ̄ - ϕᵣ) / f.time_scale
return ifelse(z < f.z_bottom, 0, tendency)
end

@glwagner
Copy link
Copy Markdown
Member

okay, then how about FieldTimeSeriesRelaxation?

Comment thread src/Forcings/fts_relaxation.jl
Comment thread src/Forcings/fts_relaxation.jl
Comment thread src/Forcings/fts_relaxation.jl Outdated
Comment on lines +199 to +201
ϕ̄ = @inbounds f.current_field[1, 1, k]
ϕᵣ = @inbounds f.target[1, 1, k]
ρ = @inbounds f.density[1, 1, k]
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.

I think this function is unnecessary. If the fields are properly reduced (Nothing locations in x, y) then they can be indexed with i, j, k.

See the getindex methods here:

https://github.com/CliMA/Oceananigans.jl/blob/19a44174e6da16e5e12d47ef56a9f9c9dc9d6a3f/src/Fields/field.jl#L570

specific_name = strip_density_prefix(name)
specific_field = context.specific_fields[specific_name]

reference_column = resolve_reference_column(forcing.reference, forcing.reference_column)
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.

@ewquon it looks like resolve_reference_column appears here

Comment thread src/Forcings/relaxation_forcing.jl Outdated
Comment thread examples/single_column_nudging.jl Outdated
Comment thread examples/single_column_nudging.jl Outdated
Comment thread examples/single_column_nudging.jl
Comment thread examples/single_column_nudging.jl Outdated
@glwagner glwagner changed the title Add RelaxationForcing Add FieldTimeSeriesRelaxation Mar 18, 2026
``z_b`` is the cutoff height below which no nudging is applied.

The `reference` [`FieldTimeSeries`](@ref) must contain the specific variable
(e.g. ``u``, ``v``, ``θ``, ``q^t``), not the density-weighted form. The user is
Copy link
Copy Markdown
Member

@giordano giordano Mar 19, 2026

Choose a reason for hiding this comment

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

Can't use @ref with symbols defined in other packages

Suggested change
(e.g. ``u``, ``v``, ``θ``, ``q^t``), not the density-weighted form. The user is
toward a reference `FieldTimeSeries`:
```math
F_{ρ ϕ} = -\\frac{ρ_r \\left(\\bar{ϕ} - ϕ_r\\right)}{τ} \\, \\mathbf{1}_{z \\geq z_b}
```
where ``\\bar{ϕ}`` is either the local field value (3D mode) or its horizontal
average (profile mode), ``ϕ_r`` is the time-interpolated reference value,
``ρ_r`` is the reference density, ``τ`` is the relaxation time scale, and
``z_b`` is the cutoff height below which no nudging is applied.
The `reference` `FieldTimeSeries` must contain the specific variable

Maybe one day we'll explore DocumenterInterLinks.jl, but we aren't using anything like that at the moment.

Edit: I think GitHub messed up my suggestion, this was meant to be [`FieldTimeSeries`](@ref) -> `FieldTimeSeries`.

ewquon and others added 10 commits March 19, 2026 23:09
Oceananigans wraps Julia functions in ContinuousBoundaryFunction{Nothing,Nothing,Nothing}
at BC construction time (before regularization). When energy_to_theta_bcs converts ρe
BCs to ρθ, it captures this unregularized ContinuousBoundaryFunction inside
EnergyFluxBoundaryConditionFunction. The second regularization pass hits the fallback
(condition → condition), leaving the inner ContinuousBoundaryFunction with all-Nothing
location parameters. This made getbc ambiguous between x-, y-, and z-boundary overloads.

Fix: add regularize_boundary_condition overloads for EnergyFluxBoundaryConditionFunction
and ThetaFluxBoundaryConditionFunction that forward regularization to the inner condition,
so the ContinuousBoundaryFunction gets proper {LX,LY,LZ} type parameters.

Add regression test for function-type ρe surface flux BCs.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
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.

I think this example should eventually be moved to NumericalEarth (and added to the docs build if we keep it). We can keep this for a future PR though.

The nudging tendency now ramps smoothly from 0 at z_bottom to 1 at z_top
via (1 - cos(π r))/2 with r = clamp((z - z_b)/(z_t - z_b), 0, 1), avoiding
the discontinuity that a hard step introduces at the boundary of the
nudged region. z_top defaults to z_bottom + 1000 m so existing call sites
keep running without changes.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
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.

4 participants