NestedSimulation abstraction + ERA5 regional hindcast example

.gitignore

@@ -68,3 +68,6 @@ CondaPkg.toml
 # claude
 .claude
 test/Manifest.toml
+
+# Throwaway scratch scripts
+scratch/

docs/make.jl

@@ -37,7 +37,8 @@ examples = [
     Example("Global climate simulation", "global_climate_simulation", false),
     # Example("Veros ocean simulation", "veros_ocean_forced_simulation", false),
     Example("Breeze over two oceans", "breeze_over_two_oceans", false),
-    Example("ERA5 hourly data", "ERA5_hourly_data", true),
+    Example("ERA5 hourly data retrieval", "ERA5_hourly_data", true),
+    Example("ERA5 downscaling with Breeze and NestedSimulation", "era5_breeze", true),
 ]
 
 # Developer examples from docs/src/developers/ directory

ext/NumericalEarthBreezeExt/NumericalEarthBreezeExt.jl

@@ -12,5 +12,6 @@ using NumericalEarth.EarthSystemModels.InterfaceComputations: ComponentExchanger
 
 include("breeze_atmosphere_interface.jl")
 include("breeze_atmosphere_simulation.jl")
+include("breeze_child_simulation.jl")
 
 end # module NumericalEarthBreezeExt

ext/NumericalEarthBreezeExt/breeze_child_simulation.jl

@@ -0,0 +1,32 @@
+#####
+##### child_simulation dispatch for Breeze.AtmosphereModel
+#####
+#
+# Two extensions to NumericalEarth.NestedSimulations:
+#
+# 1. `parent_variables(::Type{<:AtmosphereModel}, ::PrescribedAtmosphere)` —
+#    declares the (child, parent) pair mapping. Returns a NamedTuple naming
+#    the parent FTSs that drive the child's momentum BCs.
+#
+#    Density convention (important): Breeze's prognostic is density-weighted
+#    momentum (ρu, ρv, ρw). The BC value Oceananigans writes into ρu at the
+#    boundary is exactly the value of the source we hand it; so the parent's
+#    velocity FTS slots are interpreted as **momentum** values (ρ̄·u). The user
+#    populating a `PrescribedAtmosphere` for use with Breeze is expected to
+#    pre-multiply by ρ̄ when calling `set!`. See [[nested-simulations-design]]
+#    and the deferred PrescribedAtmosphere refactor (issue #266) for the
+#    slot-naming followup.
+#
+# 2. `_build_child_model(::Type{<:AtmosphereModel}, grid; …)` — routes
+#    `child_simulation(AtmosphereModel, …)` through `atmosphere_simulation`
+#    instead of constructing `AtmosphereModel` directly, so Breeze's default
+#    advection / microphysics settings apply unless the user overrides them.
+
+import NumericalEarth.EarthSystemModels.NestedSimulations: parent_variables,
+                                                            _build_child_model
+
+parent_variables(::Type{<:Breeze.AtmosphereModel}, parent::NumericalEarth.PrescribedAtmosphere) =
+    (ρu = parent.velocities.u, ρv = parent.velocities.v)
+
+_build_child_model(::Type{<:Breeze.AtmosphereModel}, grid; kwargs...) =
+    NumericalEarth.Atmospheres.atmosphere_simulation(grid; kwargs...)

src/Atmospheres/prescribed_atmosphere.jl

@@ -29,17 +29,36 @@ function Base.show(io::IO, pa::PrescribedAtmosphere)
     print(io, "└── boundary_layer_height: ", prettysummary(pa.boundary_layer_height))
 end
 
+# A grid is treated as volumetric (3D) when it carries more than one vertical
+# cell. Single-layer grids (Nz=1) are the conventional "surface atmosphere" used
+# for coupling to oceans; multi-layer grids drive the volumetric defaults used
+# when a `PrescribedAtmosphere` plays the role of a nesting parent.
+@inline is_volumetric_atmosphere_grid(grid) = size(grid, 3) > 1
+
 function default_atmosphere_velocities(grid, times)
-    ua = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    va = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    return (u=ua, v=va)
+    if is_volumetric_atmosphere_grid(grid)
+        ua = FieldTimeSeries{Center, Center, Center}(grid, times)
+        va = FieldTimeSeries{Center, Center, Center}(grid, times)
+        wa = FieldTimeSeries{Center, Center, Center}(grid, times)
+        return (u=ua, v=va, w=wa)
+    else
+        ua = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+        va = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+        return (u=ua, v=va)
+    end
 end
 
 function default_atmosphere_tracers(grid, times)
-    Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    qa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    parent(Ta) .= 273.15 + 20
-    return (T=Ta, q=qa)
+    if is_volumetric_atmosphere_grid(grid)
+        Ta = FieldTimeSeries{Center, Center, Center}(grid, times)
+        qa = FieldTimeSeries{Center, Center, Center}(grid, times)
+        return (T=Ta, q=qa)
+    else
+        Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+        qa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+        parent(Ta) .= 273.15 + 20
+        return (T=Ta, q=qa)
+    end
 end
 
 """
@@ -64,7 +83,11 @@ PrescribedPrecipitationFlux(; rain=nothing, snow=nothing) =
 Adapt.adapt_structure(to, ff::PrescribedPrecipitationFlux) =
     PrescribedPrecipitationFlux(adapt(to, ff.rain), adapt(to, ff.snow))
 
+# Surface freshwater fluxes are meaningless for a volumetric atmosphere acting as
+# a nesting parent; default to `nothing` there. The surface_*_flux accessors and
+# `extract_field_time_series` already handle the Nothing branch.
 function default_freshwater_flux(grid, times)
+    is_volumetric_atmosphere_grid(grid) && return nothing
     rain = FieldTimeSeries{Center, Center, Nothing}(grid, times)
     snow = FieldTimeSeries{Center, Center, Nothing}(grid, times)
     return PrescribedPrecipitationFlux(rain, snow)
@@ -84,9 +107,13 @@ end
 """ The standard unit of atmospheric pressure; 1 standard atmosphere (atm) = 101,325 Pascals (Pa)
 in SI units. This is approximately equal to the mean sea-level atmospheric pressure on Earth. """
 function default_atmosphere_pressure(grid, times)
-    pa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    parent(pa) .= 101325
-    return pa
+    if is_volumetric_atmosphere_grid(grid)
+        return FieldTimeSeries{Center, Center, Center}(grid, times)
+    else
+        pa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+        parent(pa) .= 101325
+        return pa
+    end
 end
 
 @inline function Oceananigans.TimeSteppers.update_state!(atmos::PrescribedAtmosphere)

src/EarthSystemModels/EarthSystemModels.jl

@@ -51,6 +51,10 @@ using .InterfaceComputations
 include("earth_system_model.jl")
 include("time_step_earth_system_model.jl")
 
+include("NestedSimulations/NestedSimulations.jl")
+
+using .NestedSimulations
+
 #####
 #####  Fallbacks for no-interface models
 #####

src/EarthSystemModels/NestedSimulations/NestedSimulations.jl

@@ -0,0 +1,24 @@
+module NestedSimulations
+
+export NestedModel, NestedSimulation,
+       parent_boundary_conditions, parent_forcings,
+       child_simulation, parent_variables
+
+using Oceananigans
+using Oceananigans.Fields: Field, Center, Face, instantiated_location
+using Oceananigans.Forcings: Relaxation
+using Oceananigans.Grids: AbstractGrid, xnode, ynode, znode
+using Oceananigans.BoundaryConditions: OpenBoundaryCondition, FieldBoundaryConditions
+using Oceananigans.OutputReaders: FieldTimeSeries
+using Oceananigans.Simulations: Simulation
+using Oceananigans.TimeSteppers: time_step!
+using Oceananigans.Units: Time
+
+include("nested_model.jl")
+include("nested_simulation.jl")
+include("interpolated_fts_boundary.jl")
+include("parent_boundary_conditions.jl")
+include("parent_forcings.jl")
+include("child_simulation.jl")
+
+end # module

src/EarthSystemModels/NestedSimulations/child_simulation.jl

@@ -0,0 +1,124 @@
+#####
+##### child_simulation: one-call construction of a parent-driven child model
+#####
+#
+# Despite the name (chosen for consistency with `ocean_simulation` /
+# `atmosphere_simulation`), `child_simulation` returns an `AbstractModel` — the
+# Breeze `atmosphere_simulation` helper returns a model, not a `Simulation`,
+# and we follow that convention so the call sites compose cleanly:
+#
+#     child  = child_simulation(modeltype, grid, parent; …)
+#     nested = NestedSimulation(parent, child; Δt, stop_time, …)
+#
+# Variable mapping (which parent FTS drives which child boundary) dispatches on
+# `parent_variables(modeltype, parent)` — multi-dispatch on the (child, parent)
+# pair. Methods live wherever both types are visible (often a package
+# extension). Cross-realm pairs (e.g. ocean child × atmosphere parent) are
+# intentionally undefined; the generic fallback throws an explanatory
+# `ArgumentError` and the user must pass `variables=` explicitly to override.
+#
+# The underlying model constructor is selected by `_build_child_model(modeltype, grid; …)`
+# — also extended via package extensions for new model types.
+
+"""
+    parent_variables(child_modeltype, parent) → NamedTuple
+
+Return the `child_variable_name => interpolatable_source` mapping naming
+the parent fields that drive each child boundary, where the source may be
+a `FieldTimeSeries` (prescribed parents) or an `AbstractField` (prognostic
+parents).
+
+Dispatch on the **(child_modeltype, parent_type)** pair — define a method
+per pair, in whichever package can see both types. The generic fallback
+errors so cross-realm pairs (e.g. an ocean child with an atmosphere parent)
+fail loudly at construction time. Users can always sidestep by passing
+`variables=` directly to [`child_simulation`](@ref).
+
+The exact contents of the NamedTuple — what fields, whether they are
+density-weighted, etc. — are determined by what the child consumes and what
+the parent provides. For example, Breeze `AtmosphereModel` consumes
+prognostic momentum (`ρu`, `ρv`, …), so a `PrescribedAtmosphere` parent
+populated for that use case must store ρ·u in its velocity slots; the
+method body documents the convention.
+"""
+function parent_variables end
+
+# Generic fallback — errors with a message that explains the situation and
+# suggests the two ways out.
+function parent_variables(child_modeltype, parent)
+    throw(ArgumentError(string(
+        "No `parent_variables` mapping is defined for ",
+        "child=$(child_modeltype) with parent=$(typeof(parent)).\n",
+        "\n",
+        "Either:\n",
+        "  • define `parent_variables(::Type{<:$(child_modeltype)}, ",
+        "parent::$(typeof(parent)))` in the package where both types are visible, or\n",
+        "  • pass `variables = (...)` explicitly to `child_simulation`.\n",
+        "\n",
+        "Note: cross-realm nesting (e.g. an ocean child driven by an atmosphere ",
+        "parent) is intentionally unsupported as a default — the variable units ",
+        "do not match. Define a custom mapping only if you own the semantics."
+    )))
+end
+
+"""
+    _build_child_model(modeltype, grid; kwargs...)
+
+Construct the child model. Defaults to `modeltype(grid; kwargs...)`; package
+extensions override this for specific model types (e.g. the Breeze ext routes
+`AtmosphereModel` through `atmosphere_simulation` so Breeze's default
+advection / microphysics apply).
+"""
+_build_child_model(modeltype, grid; kwargs...) = modeltype(grid; kwargs...)
+
+"""
+    child_simulation(modeltype, grid, parent;
+                     sides = (:west, :east, :south, :north),
+                     schemes = NamedTuple(),
+                     variables = parent_variables(modeltype, parent),
+                     relaxation_rate = nothing,
+                     relaxation_mask = 1,
+                     model_kwargs...)
+
+Construct a child model (of type `modeltype`) on `grid` with Open boundary
+conditions driven by `parent`. Returns the constructed `AbstractModel`.
+
+`parent_boundary_conditions` wires the BCs from `variables`; `model_kwargs`
+are forwarded to the underlying constructor (see `_build_child_model`).
+
+If `relaxation_rate` is supplied, [`parent_forcings`](@ref) builds an
+Oceananigans `Relaxation` per variable that nudges the child interior toward
+the parent at rate `relaxation_rate`, weighted by `relaxation_mask` (default
+uniform). Per-variable values are supported via `NamedTuple`s. The resulting
+forcings are merged with any `forcing` already passed in `model_kwargs`.
+
+Wrap with `NestedSimulation(parent, child; Δt, stop_time, …)` to integrate.
+"""
+function child_simulation(modeltype, grid, parent;
+                          sides = (:west, :east, :south, :north),
+                          schemes = NamedTuple(),
+                          variables = parent_variables(modeltype, parent),
+                          relaxation_rate = nothing,
+                          relaxation_mask = 1,
+                          model_kwargs...)
+
+    bcs = parent_boundary_conditions(grid; variables, sides, schemes)
+
+    # Pull any user-supplied `forcing` out so we can merge in the relaxation
+    # forcings on the same field names without colliding.
+    model_kwargs_nt = NamedTuple(model_kwargs)
+    user_forcing   = get(model_kwargs_nt, :forcing, NamedTuple())
+    rest_kwargs    = Base.structdiff(model_kwargs_nt, NamedTuple{(:forcing,)})
+
+    forcing = if relaxation_rate === nothing
+        user_forcing
+    else
+        relaxation = parent_forcings(; variables, rate = relaxation_rate, mask = relaxation_mask)
+        merge(user_forcing, relaxation)
+    end
+
+    return _build_child_model(modeltype, grid;
+                              boundary_conditions = bcs,
+                              forcing,
+                              rest_kwargs...)
+end