From cd250d9d65066a31cb9c2e48a1f1680004fffba2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:07:47 +0200 Subject: [PATCH 01/13] Replace XESMF with ConservativeRegridding as hard dependency Remove XESMF entirely and use ConservativeRegridding for regridding between SpeedyWeather and Oceananigans grids. ConservativeRegridding is now a direct dependency rather than an extension. The SpeedyWeather extension now implements treeify for RingGrids, building GeoInterface polygons from cell vertices. This enables the default OctahedralGaussianGrid instead of requiring FullClenshawGrid. Note: treeify for SpeedyWeather grids needs testing to confirm it works end-to-end with the conservative regridding pipeline. Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 6 +- docs/Project.toml | 1 - examples/global_climate_simulation.jl | 6 +- .../NumericalEarthSpeedyWeatherExt.jl | 6 +- .../speedy_regridder.jl | 76 +++++++------------ .../speedy_weather_exchanger.jl | 71 +++++++++-------- test/Project.toml | 2 - test/runtests_setup.jl | 8 -- test/test_speedy_coupling.jl | 4 +- 9 files changed, 77 insertions(+), 103 deletions(-) diff --git a/Project.toml b/Project.toml index 0f01f8982..7ce7d5615 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ authors = ["NumericalEarth contributors"] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CFTime = "179af706-886a-5703-950a-314cd64e0468" ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" +ConservativeRegridding = "8e50ac2c-eb48-49bc-a402-07c87b949343" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -39,14 +40,13 @@ PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Reactant = "3c362404-f566-11ee-1572-e11a4b42c853" SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9" WorldOceanAtlasTools = "04f20302-f1b9-11e8-29d9-7d841cb0a64a" -XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" [extensions] NumericalEarthBreezeExt = "Breeze" NumericalEarthCDSAPIExt = "CDSAPI" NumericalEarthCopernicusMarineExt = "CopernicusMarine" NumericalEarthReactantExt = "Reactant" -NumericalEarthSpeedyWeatherExt = ["SpeedyWeather", "XESMF"] +NumericalEarthSpeedyWeatherExt = "SpeedyWeather" NumericalEarthVerosExt = ["PythonCall", "CondaPkg"] NumericalEarthWOAExt = "WorldOceanAtlasTools" @@ -80,7 +80,7 @@ SpeedyWeather = "0.17.4" StaticArrays = "1" Statistics = "<0.0.1, 1" Thermodynamics = "0.15.3" +ConservativeRegridding = "0.2" WorldOceanAtlasTools = "0.6" -XESMF = "0.1.6" ZipFile = "0.10" julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml index c10ac23e9..f0eab7673 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -17,7 +17,6 @@ Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9" -XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" [sources] NumericalEarth = {path = ".."} diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl index 5a40eeac2..b0fa5fbea 100644 --- a/examples/global_climate_simulation.jl +++ b/examples/global_climate_simulation.jl @@ -17,9 +17,9 @@ # `SpeedyWeather.PrimitiveWetModel` (the atmosphere). All these are coupled and orchestrated by the `NumericalEarth.EarthSystemModel` # (the coupled system). # -# The XESMF.jl package is used to regrid fields between the atmosphere and ocean--sea ice components. +# The ConservativeRegridding.jl package is used to regrid fields between the atmosphere and ocean--sea ice components. -using Oceananigans, SpeedyWeather, XESMF, NumericalEarth +using Oceananigans, SpeedyWeather, NumericalEarth using NCDatasets, CairoMakie using Oceananigans.Units using Printf, Statistics, Dates @@ -76,7 +76,7 @@ Oceananigans.set!(sea_ice.model, h=Metadatum(:sea_ice_thickness, dataset=ECCO4Mo # The `atmosphere_simulation` function takes care of building an atmosphere model with appropriate # hooks so that NumericalEarth can compute inter-component fluxes. nlayers = 4 -spectral_grid = SpeedyWeather.SpectralGrid(; trunc=63, nlayers, Grid=FullClenshawGrid) +spectral_grid = SpeedyWeather.SpectralGrid(; trunc=63, nlayers) atmosphere = atmosphere_simulation(spectral_grid, output=true) # The atmosphere model already includes some initial conditions as described above: diff --git a/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl b/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl index a583861fb..047295a94 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl @@ -4,9 +4,9 @@ using OffsetArrays using KernelAbstractions using Statistics -import SpeedyWeather -import NumericalEarth -import Oceananigans +import SpeedyWeather +import NumericalEarth +import Oceananigans import SpeedyWeather.RingGrids include("speedy_atmosphere_simulations.jl") diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl index 85b2db32c..a1c44274e 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl @@ -1,49 +1,29 @@ -using Oceananigans.Grids: AbstractGrid -using Oceananigans - -import XESMF: Regridder, xesmf_coordinates - -const Grids = Union{SpeedyWeather.SpectralGrid, AbstractGrid} - -function Regridder(src::Grids, dst::Grids; method::String="bilinear", periodic=true) - src_coords = xesmf_coordinates(src, Center(), Center(), Center()) - dst_coords = xesmf_coordinates(dst, Center(), Center(), Center()) - - return XESMF.Regridder(src_coords, dst_coords; method, periodic) -end - -two_dimensionalize(lat::Matrix, lon::Matrix) = lat, lon - -function two_dimensionalize(lat::AbstractVector, lon::AbstractVector) - Nx = length(lon) - Ny = length(lat) - lat = repeat(lat', Nx) - lon = repeat(lon, 1, Ny) - return lat, lon -end - -xesmf_coordinates(grid::SpeedyWeather.SpectralGrid, args...) = xesmf_coordinates(grid.grid, args...) -xesmf_coordinates(grid::SpeedyWeather.RingGrids.AbstractGrid, args...) = - throw(ArgumentError("xesmf_coordinates not implemented for grid type $(typeof(grid)), maybe you meant to pass a FullGrid?")) - -function xesmf_coordinates(grid::SpeedyWeather.RingGrids.AbstractFullGrid, args...) - lon = RingGrids.get_lond(grid) - lat = RingGrids.get_latd(grid) - dlon = lon[2] - lon[1] - - lat_b = [90, 0.5 .* (lat[1:end-1] .+ lat[2:end])..., -90] - lon_b = [lon[1] - dlon / 2, lon .+ dlon / 2...] - - lat, lon = two_dimensionalize(lat, lon) - lat_b, lon_b = two_dimensionalize(lat_b, lon_b) - - # Python's xESMF expects 2D arrays with (x, y) coordinates - # in which y varies in dim=1 and x varies in dim=2 - # therefore we transpose the coordinate matrices - coords_dictionary = Dict("lat" => permutedims(lat, (2, 1)), # φ is latitude - "lon" => permutedims(lon, (2, 1)), # λ is longitude - "lat_b" => permutedims(lat_b, (2, 1)), - "lon_b" => permutedims(lon_b, (2, 1))) - - return coords_dictionary +using ConservativeRegridding: Trees +import ConservativeRegridding.Trees: treeify +import GeometryOpsCore as GOCore +import GeometryOps as GO +import GeoInterface as GI + +using StaticArrays: SA + +GOCore.best_manifold(grid::RingGrids.AbstractGrid) = GO.Spherical() +GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.grid) + +treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(manifold, sg.grid) + +function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) + polygons_matrix = RingGrids.get_gridcell_polygons(grid) + npoints = size(polygons_matrix, 2) + lonlat_to_usp = GO.UnitSphereFromGeographic() + + polys = Vector{GI.Polygon}(undef, npoints) + for ij in 1:npoints + v1 = lonlat_to_usp(polygons_matrix[1, ij]) + v2 = lonlat_to_usp(polygons_matrix[2, ij]) + v3 = lonlat_to_usp(polygons_matrix[3, ij]) + v4 = lonlat_to_usp(polygons_matrix[4, ij]) + polys[ij] = GI.Polygon(SA[GI.LinearRing(SA[v1, v2, v3, v4, v1])]) + end + + return Trees.treeify(manifold, polys) end diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl index 08c024e71..cc2b86ed7 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl @@ -3,13 +3,14 @@ using Oceananigans.BoundaryConditions using Oceananigans.Grids: architecture using Oceananigans.Utils: launch! using Oceananigans.Operators: intrinsic_vector -using XESMF + +using ConservativeRegridding: Regridder +using LinearAlgebra: transpose + +import ConservativeRegridding: regrid! using NumericalEarth.EarthSystemModels: sea_ice_concentration -# TODO: Implement conservative regridding when ready -# using ConservativeRegridding -# using GeoInterface: Polygon, LinearRing import NumericalEarth.EarthSystemModels: update_net_fluxes!, interpolate_state! import NumericalEarth.Atmospheres: atmosphere_regridder import NumericalEarth.EarthSystemModels.InterfaceComputations: net_fluxes, ComponentExchanger @@ -23,14 +24,13 @@ net_fluxes(::SpeedySimulation) = nothing # If this work we can # 1. Copy speedyweather gridarrays to the GPU # 2. Perform the regridding on the GPU -function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) +function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) spectral_grid = atmosphere.model.spectral_grid - # TODO: Implement a conservative regridder when ready - from_atmosphere = XESMF.Regridder(spectral_grid, exchange_grid) - to_atmosphere = XESMF.Regridder(exchange_grid, spectral_grid) + from_atmosphere = Regridder(exchange_grid, spectral_grid) + to_atmosphere = transpose(from_atmosphere) regridder = (; to_atmosphere, from_atmosphere) - + state = (; u = Field{Center, Center, Nothing}(exchange_grid), v = Field{Center, Center, Nothing}(exchange_grid), T = Field{Center, Center, Nothing}(exchange_grid), @@ -39,18 +39,23 @@ function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) ℐꜜˢʷ = Field{Center, Center, Nothing}(exchange_grid), ℐꜜˡʷ = Field{Center, Center, Nothing}(exchange_grid), Jᶜ = Field{Center, Center, Nothing}(exchange_grid)) - + return ComponentExchanger(state, regridder) end -@inline (regrid!::XESMF.Regridder)(field::Oceananigans.Field, data::AbstractArray) = regrid!(vec(interior(field)), data) -@inline (regrid!::XESMF.Regridder)(data::AbstractArray, field::Oceananigans.Field) = regrid!(data, vec(interior(field))) +function regrid!(field::Oceananigans.Field, regridder::Regridder, data::AbstractArray) + regrid!(vec(interior(field)), regridder, vec(data)) +end + +function regrid!(data::AbstractArray, regridder::Regridder, field::Oceananigans.Field) + regrid!(vec(data), regridder, vec(interior(field))) +end # Regrid the atmospheric state on the exchange grid function interpolate_state!(exchanger, exchange_grid, atmos::SpeedySimulation, coupled_model) - regrid! = exchanger.regridder.from_atmosphere - exchange_state = exchanger.state - surface_layer = atmos.model.spectral_grid.nlayers + from_atmosphere = exchanger.regridder.from_atmosphere + exchange_state = exchanger.state + surface_layer = atmos.model.spectral_grid.nlayers ua = RingGrids.field_view(atmos.diagnostic_variables.grid.u_grid, :, surface_layer).data va = RingGrids.field_view(atmos.diagnostic_variables.grid.v_grid, :, surface_layer).data @@ -61,14 +66,14 @@ function interpolate_state!(exchanger, exchange_grid, atmos::SpeedySimulation, c ℐꜜˡʷ = atmos.diagnostic_variables.physics.surface_longwave_down.data Jᶜ = atmos.diagnostic_variables.physics.total_precipitation_rate.data - regrid!(exchange_state.u, ua) - regrid!(exchange_state.v, va) - regrid!(exchange_state.T, Ta) - regrid!(exchange_state.q, qa) - regrid!(exchange_state.p, pa) - regrid!(exchange_state.ℐꜜˢʷ, ℐꜜˢʷ) - regrid!(exchange_state.ℐꜜˡʷ, ℐꜜˡʷ) - regrid!(exchange_state.Jᶜ, Jᶜ) + regrid!(exchange_state.u, from_atmosphere, ua) + regrid!(exchange_state.v, from_atmosphere, va) + regrid!(exchange_state.T, from_atmosphere, Ta) + regrid!(exchange_state.q, from_atmosphere, qa) + regrid!(exchange_state.p, from_atmosphere, pa) + regrid!(exchange_state.ℐꜜˢʷ, from_atmosphere, ℐꜜˢʷ) + regrid!(exchange_state.ℐꜜˡʷ, from_atmosphere, ℐꜜˡʷ) + regrid!(exchange_state.Jᶜ, from_atmosphere, Jᶜ) arch = architecture(exchange_grid) @@ -90,17 +95,17 @@ end @kernel function _rotate_winds!(u, v, grid) i, j = @index(Global, NTuple) - kᴺ = size(grid, 3) + kᴺ = size(grid, 3) uₑ, vₑ = intrinsic_vector(i, j, kᴺ, grid, u, v) @inbounds u[i, j, kᴺ] = uₑ @inbounds v[i, j, kᴺ] = vₑ end -# TODO: Fix the coupling with the sea ice model and make sure that +# TODO: Fix the coupling with the sea ice model and make sure that # the this function works also for sea_ice=nothing and on GPUs without # needing to allocate memory. function update_net_fluxes!(coupled_model, atmos::SpeedySimulation) - regrid! = coupled_model.interfaces.exchanger.atmosphere.regridder.to_atmosphere + to_atmosphere = coupled_model.interfaces.exchanger.atmosphere.regridder.to_atmosphere ao_fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes ai_fluxes = coupled_model.interfaces.atmosphere_sea_ice_interface.fluxes @@ -120,16 +125,16 @@ function update_net_fluxes!(coupled_model, atmos::SpeedySimulation) # TODO: Figure out how we are going to deal with upwelling radiation # TODO: regrid longwave rather than a mixed surface temperature # TODO: This does not work on GPUs!! - regrid!(𝒬ᵀ_speedy, vec(interior(𝒬ᵀᵃᵒ) .* (1 .- ℵ) .+ ℵ .* interior(𝒬ᵀᵃⁱ))) - regrid!(Jᵛ_speedy, vec(interior(Jᵛᵃᵒ) .* (1 .- ℵ) .+ ℵ .* interior(Jᵛᵃⁱ))) - regrid!(sst, vec(interior(To) .* (1 .- ℵ) .+ ℵ .* interior(Ti) .+ 273.15)) + regrid!(𝒬ᵀ_speedy, to_atmosphere, vec(interior(𝒬ᵀᵃᵒ) .* (1 .- ℵ) .+ ℵ .* interior(𝒬ᵀᵃⁱ))) + regrid!(Jᵛ_speedy, to_atmosphere, vec(interior(Jᵛᵃᵒ) .* (1 .- ℵ) .+ ℵ .* interior(Jᵛᵃⁱ))) + regrid!(sst, to_atmosphere, vec(interior(To) .* (1 .- ℵ) .+ ℵ .* interior(Ti) .+ 273.15)) return nothing end # Simple case -> there is no sea ice! function update_net_fluxes!(coupled_model::SpeedyNoSeaIceEarthSystemModel, atmos::SpeedySimulation) - regrid! = coupled_model.interfaces.exchanger.atmosphere.regridder.to_atmosphere + to_atmosphere = coupled_model.interfaces.exchanger.atmosphere.regridder.to_atmosphere ao_fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes 𝒬ᵀᵃᵒ = ao_fluxes.sensible_heat Jᵛᵃᵒ = ao_fluxes.water_vapor @@ -142,9 +147,9 @@ function update_net_fluxes!(coupled_model::SpeedyNoSeaIceEarthSystemModel, atmos # TODO: Figure out how we are going to deal with upwelling radiation # TODO: This does not work on GPUs!! - regrid!(𝒬ᵀ_speedy, vec(interior(𝒬ᵀᵃᵒ))) - regrid!(Jᵛ_speedy, vec(interior(Jᵛᵃᵒ))) - regrid!(sst, vec(interior(To) .+ 273.15)) + regrid!(𝒬ᵀ_speedy, to_atmosphere, vec(interior(𝒬ᵀᵃᵒ))) + regrid!(Jᵛ_speedy, to_atmosphere, vec(interior(Jᵛᵃᵒ))) + regrid!(sst, to_atmosphere, vec(interior(To) .+ 273.15)) return nothing end diff --git a/test/Project.toml b/test/Project.toml index 8b9e4d34e..7b69366c2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -25,7 +25,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" WorldOceanAtlasTools = "04f20302-f1b9-11e8-29d9-7d841cb0a64a" -XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" [sources] NumericalEarth = {path = ".."} @@ -55,7 +54,6 @@ Statistics = "<0.0.1, 1" Test = "<0.0.1, 1" Thermodynamics = "0.15.3" WorldOceanAtlasTools = "0.6" -XESMF = "0.1.6" julia = "1.10" [extras] diff --git a/test/runtests_setup.jl b/test/runtests_setup.jl index b3e0fae9d..03909ff32 100644 --- a/test/runtests_setup.jl +++ b/test/runtests_setup.jl @@ -50,14 +50,6 @@ test_fields = Dict( EN4Monthly() => (:T, :S), ) -# Install XESMF for windows architectures -if Sys.iswindows() && Sys.ARCH == :x86_64 - @info "Installing XESMF for windows architectures" - using Pkg - Pkg.add("CondaPkg") - using CondaPkg - CondaPkg.add(["esmf", "esmpy"]) -end ##### ##### Test utilities diff --git a/test/test_speedy_coupling.jl b/test/test_speedy_coupling.jl index 5a234001c..772141589 100644 --- a/test/test_speedy_coupling.jl +++ b/test/test_speedy_coupling.jl @@ -1,4 +1,4 @@ -using SpeedyWeather, XESMF +using SpeedyWeather using NumericalEarth using Oceananigans using Dates @@ -7,7 +7,7 @@ using Test NumericalEarthSpeedyWeatherExt = Base.get_extension(NumericalEarth, :NumericalEarthSpeedyWeatherExt) @test !isnothing(NumericalEarthSpeedyWeatherExt) -spectral_grid = SpeedyWeather.SpectralGrid(trunc=51, nlayers=3, Grid=FullClenshawGrid) +spectral_grid = SpeedyWeather.SpectralGrid(trunc=51, nlayers=3) oceananigans_grid = LatitudeLongitudeGrid(Oceananigans.CPU(); size=(200, 100, 1), latitude=(-80, 80), longitude=(0, 360), z = (0, 1)) ocean = NumericalEarth.Oceans.ocean_simulation(oceananigans_grid; momentum_advection=nothing, tracer_advection=nothing, closure=nothing) From 5b55493870069817cb2679065e7111a4d566a368 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:10:20 +0200 Subject: [PATCH 02/13] run example in PR --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index f3f332f8d..66199adf3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -32,7 +32,7 @@ examples = [ Example("Single-column ocean simulation", "single_column_os_papa_simulation", true), Example("One-degree ocean--sea ice simulation", "one_degree_simulation", false), Example("Near-global ocean simulation", "near_global_ocean_simulation", false), - Example("Global climate simulation", "global_climate_simulation", false), + Example("Global climate simulation", "global_climate_simulation", true), Example("Veros ocean simulation", "veros_ocean_forced_simulation", false), Example("Breeze over two oceans", "breeze_over_two_oceans", false), Example("ERA5 winds and Stokes drift", "ERA5_winds_and_stokes_drift", true), From 03b5ce4f413baf7a046de1889b2fe6650114a7e5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:10:41 +0200 Subject: [PATCH 03/13] Fix vertex ordering in treeify to counter-clockwise RingGrids.get_gridcell_polygons returns vertices in clockwise order (E, S, W, N). Reverse to CCW (E, N, W, S) so that GO.area on the Spherical manifold returns positive areas. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../speedy_regridder.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl index a1c44274e..209333328 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl @@ -12,17 +12,21 @@ GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.g treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(manifold, sg.grid) function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) + # get_gridcell_polygons returns a 5×npoints matrix of (lon, lat) tuples + # with vertices in clockwise order: E, S, W, N, E (closed). + # We reverse to counter-clockwise (E, N, W, S, E) for correct signed area. polygons_matrix = RingGrids.get_gridcell_polygons(grid) npoints = size(polygons_matrix, 2) lonlat_to_usp = GO.UnitSphereFromGeographic() polys = Vector{GI.Polygon}(undef, npoints) for ij in 1:npoints - v1 = lonlat_to_usp(polygons_matrix[1, ij]) - v2 = lonlat_to_usp(polygons_matrix[2, ij]) - v3 = lonlat_to_usp(polygons_matrix[3, ij]) - v4 = lonlat_to_usp(polygons_matrix[4, ij]) - polys[ij] = GI.Polygon(SA[GI.LinearRing(SA[v1, v2, v3, v4, v1])]) + vE = lonlat_to_usp(polygons_matrix[1, ij]) + vS = lonlat_to_usp(polygons_matrix[2, ij]) + vW = lonlat_to_usp(polygons_matrix[3, ij]) + vN = lonlat_to_usp(polygons_matrix[4, ij]) + # CCW: E → N → W → S → E + polys[ij] = GI.Polygon(SA[GI.LinearRing(SA[vE, vN, vW, vS, vE])]) end return Trees.treeify(manifold, polys) From d98db40f60db17e8e58493960e569ca8e4884613 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:17:03 +0200 Subject: [PATCH 04/13] Add quadtree spatial acceleration for RingGrids treeify Build a proper ExplicitPolygonGrid + TopDownQuadtreeCursor for SpeedyWeather grids instead of using FlatNoTree (brute force). - Full grids: clean (nlon, nlat) ExplicitPolygonGrid, linear indices match RingGrid flat indices directly. - Reduced grids: (nlon_max, nlat) matrix with degenerate padding for shorter rings, ReorderedTopDownQuadtreeCursor maps 2D indices to RingGrid flat indices, ghost cells beyond npoints have zero area. The regrid! overloads handle the size mismatch for reduced grids by zero-padding into temp arrays. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../speedy_regridder.jl | 81 ++++++++++++++++--- .../speedy_weather_exchanger.jl | 26 +++++- 2 files changed, 92 insertions(+), 15 deletions(-) diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl index 209333328..967270f13 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl @@ -1,4 +1,4 @@ -using ConservativeRegridding: Trees +using ConservativeRegridding: Trees, Regridder import ConservativeRegridding.Trees: treeify import GeometryOpsCore as GOCore import GeometryOps as GO @@ -11,23 +11,78 @@ GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.g treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(manifold, sg.grid) +function _make_degenerate_polygon(lonlat_to_usp) + p = lonlat_to_usp((0.0, 0.0)) + return GI.Polygon(SA[GI.LinearRing(SA[p, p, p, p, p])]) +end + +function _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) + vE = lonlat_to_usp(polygons_matrix[1, ij]) + vS = lonlat_to_usp(polygons_matrix[2, ij]) + vW = lonlat_to_usp(polygons_matrix[3, ij]) + vN = lonlat_to_usp(polygons_matrix[4, ij]) + # get_gridcell_polygons returns CW (E, S, W, N); reverse to CCW (E, N, W, S) + return GI.Polygon(SA[GI.LinearRing(SA[vE, vN, vW, vS, vE])]) +end + +# Full grids: all rings have the same nlon, so no padding needed. +# Linear index in ExplicitPolygonGrid matches the RingGrid flat index. +function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractFullGrid) + polygons_matrix = RingGrids.get_gridcell_polygons(grid) + lonlat_to_usp = GO.UnitSphereFromGeographic() + + nlat = RingGrids.get_nlat(grid) + nlon = RingGrids.get_nlon(grid) + + poly_matrix = [_to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, (j - 1) * nlon + i) + for i in 1:nlon, j in 1:nlat] + + epg = Trees.ExplicitPolygonGrid(manifold, poly_matrix) + tree = Trees.TopDownQuadtreeCursor(epg) + return Trees.KnownFullSphereExtentWrapper(tree) +end + +# Reduced grids: variable nlon per ring, pad shorter rings with degenerate polygons. +# The ReorderedTopDownQuadtreeCursor remaps 2D (i,j) → flat RingGrid indices 1:npoints. +# Ghost cells (padding) get indices npoints+1:ntotal and have zero intersection area. function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) - # get_gridcell_polygons returns a 5×npoints matrix of (lon, lat) tuples - # with vertices in clockwise order: E, S, W, N, E (closed). - # We reverse to counter-clockwise (E, N, W, S, E) for correct signed area. polygons_matrix = RingGrids.get_gridcell_polygons(grid) npoints = size(polygons_matrix, 2) lonlat_to_usp = GO.UnitSphereFromGeographic() - polys = Vector{GI.Polygon}(undef, npoints) - for ij in 1:npoints - vE = lonlat_to_usp(polygons_matrix[1, ij]) - vS = lonlat_to_usp(polygons_matrix[2, ij]) - vW = lonlat_to_usp(polygons_matrix[3, ij]) - vN = lonlat_to_usp(polygons_matrix[4, ij]) - # CCW: E → N → W → S → E - polys[ij] = GI.Polygon(SA[GI.LinearRing(SA[vE, vN, vW, vS, vE])]) + rings = RingGrids.eachring(grid) + nlat = length(rings) + nlon_max = RingGrids.get_nlon_max(grid) + ntotal = nlon_max * nlat + + degenerate = _make_degenerate_polygon(lonlat_to_usp) + poly_type = typeof(degenerate) + poly_matrix = Matrix{poly_type}(undef, nlon_max, nlat) + + # Reordering: real cells → RingGrid flat index (1:npoints) + # Ghost cells → npoints+1:ntotal + cart2lin = zeros(Int, nlon_max, nlat) + lin2cart = Vector{CartesianIndex{2}}(undef, ntotal) + + n_padding = 0 + for (j, ring) in enumerate(rings) + nlon_ring = length(ring) + for (i_local, ij) in enumerate(ring) + poly_matrix[i_local, j] = _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) + cart2lin[i_local, j] = ij + lin2cart[ij] = CartesianIndex(i_local, j) + end + for i_pad in (nlon_ring + 1):nlon_max + poly_matrix[i_pad, j] = degenerate + n_padding += 1 + pad_idx = npoints + n_padding + cart2lin[i_pad, j] = pad_idx + lin2cart[pad_idx] = CartesianIndex(i_pad, j) + end end - return Trees.treeify(manifold, polys) + epg = Trees.ExplicitPolygonGrid(manifold, poly_matrix) + ordering = Trees.Reorderer2D(cart2lin, lin2cart) + tree = Trees.ReorderedTopDownQuadtreeCursor(epg, ordering) + return Trees.KnownFullSphereExtentWrapper(tree) end diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl index cc2b86ed7..0220a56b5 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl @@ -43,12 +43,34 @@ function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) return ComponentExchanger(state, regridder) end +# Regrid from atmosphere (src) onto an Oceananigans field (dst). +# For reduced grids the regridder has ghost cells: src_areas is longer +# than the actual data vector, so we zero-pad into src_temp. function regrid!(field::Oceananigans.Field, regridder::Regridder, data::AbstractArray) - regrid!(vec(interior(field)), regridder, vec(data)) + src = vec(data) + nsrc = length(regridder.src_areas) + if length(src) < nsrc + fill!(regridder.src_temp, 0) + regridder.src_temp[1:length(src)] .= src + regrid!(regridder.dst_temp, regridder, regridder.src_temp) + vec(interior(field)) .= regridder.dst_temp + else + regrid!(vec(interior(field)), regridder, src) + end end +# Regrid from an Oceananigans field (src) onto atmosphere data (dst). +# For reduced grids the regridder's dst has ghost cells, so we regrid +# into dst_temp and copy back only the real entries. function regrid!(data::AbstractArray, regridder::Regridder, field::Oceananigans.Field) - regrid!(vec(data), regridder, vec(interior(field))) + dst = vec(data) + ndst = length(regridder.dst_areas) + if length(dst) < ndst + regrid!(regridder.dst_temp, regridder, regridder.src_temp .= vec(interior(field))) + dst .= view(regridder.dst_temp, 1:length(dst)) + else + regrid!(dst, regridder, vec(interior(field))) + end end # Regrid the atmospheric state on the exchange grid From a67f3da48905264393417cb8cd84df1046e0f0f5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:21:05 +0200 Subject: [PATCH 05/13] Use MultiTreeWrapper for reduced grids, revert regrid! to simple form MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Reduced grids now use one sub-tree per latitude ring combined via MultiTreeWrapper. Each ring is an ExplicitPolygonGrid(nlon_ring, 1) with IndexOffsetQuadtreeCursor mapping to flat RingGrid indices. This gives ncells = npoints exactly — no ghost cells, no size mismatch, so regrid! stays clean with no padding logic. The dual DFS prunes entire latitude rings via SphericalCap extents, and within each ring prunes by longitude via the quadtree. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../speedy_regridder.jl | 58 +++++++------------ .../speedy_weather_exchanger.jl | 26 +-------- 2 files changed, 23 insertions(+), 61 deletions(-) diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl index 967270f13..79a29c344 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl @@ -11,11 +11,6 @@ GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.g treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(manifold, sg.grid) -function _make_degenerate_polygon(lonlat_to_usp) - p = lonlat_to_usp((0.0, 0.0)) - return GI.Polygon(SA[GI.LinearRing(SA[p, p, p, p, p])]) -end - function _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) vE = lonlat_to_usp(polygons_matrix[1, ij]) vS = lonlat_to_usp(polygons_matrix[2, ij]) @@ -42,47 +37,36 @@ function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractFullGrid) return Trees.KnownFullSphereExtentWrapper(tree) end -# Reduced grids: variable nlon per ring, pad shorter rings with degenerate polygons. -# The ReorderedTopDownQuadtreeCursor remaps 2D (i,j) → flat RingGrid indices 1:npoints. -# Ghost cells (padding) get indices npoints+1:ntotal and have zero intersection area. +# Reduced grids: variable nlon per ring. Build one sub-tree per latitude ring +# and combine with MultiTreeWrapper. Each ring is an ExplicitPolygonGrid(nlon_ring, 1) +# wrapped in IndexOffsetQuadtreeCursor so that emitted indices match RingGrid flat indices. +# ncells = npoints exactly — no ghost cells, no size mismatch in regrid!. function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) polygons_matrix = RingGrids.get_gridcell_polygons(grid) - npoints = size(polygons_matrix, 2) lonlat_to_usp = GO.UnitSphereFromGeographic() - rings = RingGrids.eachring(grid) - nlat = length(rings) - nlon_max = RingGrids.get_nlon_max(grid) - ntotal = nlon_max * nlat + rings = RingGrids.eachring(grid) + nlat = length(rings) - degenerate = _make_degenerate_polygon(lonlat_to_usp) - poly_type = typeof(degenerate) - poly_matrix = Matrix{poly_type}(undef, nlon_max, nlat) + # Build one sub-tree per latitude ring + poly_type = typeof(_to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, 1)) - # Reordering: real cells → RingGrid flat index (1:npoints) - # Ghost cells → npoints+1:ntotal - cart2lin = zeros(Int, nlon_max, nlat) - lin2cart = Vector{CartesianIndex{2}}(undef, ntotal) + subtrees = Trees.IndexOffsetQuadtreeCursor[] + offsets = Int[] + cumulative = 0 - n_padding = 0 - for (j, ring) in enumerate(rings) + for ring in rings nlon_ring = length(ring) - for (i_local, ij) in enumerate(ring) - poly_matrix[i_local, j] = _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) - cart2lin[i_local, j] = ij - lin2cart[ij] = CartesianIndex(i_local, j) - end - for i_pad in (nlon_ring + 1):nlon_max - poly_matrix[i_pad, j] = degenerate - n_padding += 1 - pad_idx = npoints + n_padding - cart2lin[i_pad, j] = pad_idx - lin2cart[pad_idx] = CartesianIndex(i_pad, j) + poly_ring = Matrix{poly_type}(undef, nlon_ring, 1) + for (i, ij) in enumerate(ring) + poly_ring[i, 1] = _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) end + epg = Trees.ExplicitPolygonGrid(manifold, poly_ring) + cursor = Trees.IndexOffsetQuadtreeCursor(epg, first(ring) - 1) + push!(subtrees, cursor) + cumulative += nlon_ring + push!(offsets, cumulative) end - epg = Trees.ExplicitPolygonGrid(manifold, poly_matrix) - ordering = Trees.Reorderer2D(cart2lin, lin2cart) - tree = Trees.ReorderedTopDownQuadtreeCursor(epg, ordering) - return Trees.KnownFullSphereExtentWrapper(tree) + return Trees.MultiTreeWrapper(subtrees, offsets) end diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl index 0220a56b5..cc2b86ed7 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl @@ -43,34 +43,12 @@ function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) return ComponentExchanger(state, regridder) end -# Regrid from atmosphere (src) onto an Oceananigans field (dst). -# For reduced grids the regridder has ghost cells: src_areas is longer -# than the actual data vector, so we zero-pad into src_temp. function regrid!(field::Oceananigans.Field, regridder::Regridder, data::AbstractArray) - src = vec(data) - nsrc = length(regridder.src_areas) - if length(src) < nsrc - fill!(regridder.src_temp, 0) - regridder.src_temp[1:length(src)] .= src - regrid!(regridder.dst_temp, regridder, regridder.src_temp) - vec(interior(field)) .= regridder.dst_temp - else - regrid!(vec(interior(field)), regridder, src) - end + regrid!(vec(interior(field)), regridder, vec(data)) end -# Regrid from an Oceananigans field (src) onto atmosphere data (dst). -# For reduced grids the regridder's dst has ghost cells, so we regrid -# into dst_temp and copy back only the real entries. function regrid!(data::AbstractArray, regridder::Regridder, field::Oceananigans.Field) - dst = vec(data) - ndst = length(regridder.dst_areas) - if length(dst) < ndst - regrid!(regridder.dst_temp, regridder, regridder.src_temp .= vec(interior(field))) - dst .= view(regridder.dst_temp, 1:length(dst)) - else - regrid!(dst, regridder, vec(interior(field))) - end + regrid!(vec(data), regridder, vec(interior(field))) end # Regrid the atmospheric state on the exchange grid From de5224e57e329d972da51bb1140fddebf0344a6d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:24:22 +0200 Subject: [PATCH 06/13] Clean up speedy_regridder: shared const, simpler helpers Co-Authored-By: Claude Opus 4.6 (1M context) --- .../speedy_regridder.jl | 40 ++++++++----------- 1 file changed, 16 insertions(+), 24 deletions(-) diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl index 79a29c344..d9180b7a3 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl @@ -6,33 +6,33 @@ import GeoInterface as GI using StaticArrays: SA +const lonlat_to_unit_sphere = GO.UnitSphereFromGeographic() + GOCore.best_manifold(grid::RingGrids.AbstractGrid) = GO.Spherical() GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.grid) treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(manifold, sg.grid) -function _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) - vE = lonlat_to_usp(polygons_matrix[1, ij]) - vS = lonlat_to_usp(polygons_matrix[2, ij]) - vW = lonlat_to_usp(polygons_matrix[3, ij]) - vN = lonlat_to_usp(polygons_matrix[4, ij]) - # get_gridcell_polygons returns CW (E, S, W, N); reverse to CCW (E, N, W, S) +# get_gridcell_polygons returns CW (E, S, W, N); reverse to CCW (E, N, W, S) +function ccw_unit_sphere_polygon(polygons_matrix, ij) + vE = lonlat_to_unit_sphere(polygons_matrix[1, ij]) + vN = lonlat_to_unit_sphere(polygons_matrix[4, ij]) + vW = lonlat_to_unit_sphere(polygons_matrix[3, ij]) + vS = lonlat_to_unit_sphere(polygons_matrix[2, ij]) return GI.Polygon(SA[GI.LinearRing(SA[vE, vN, vW, vS, vE])]) end # Full grids: all rings have the same nlon, so no padding needed. # Linear index in ExplicitPolygonGrid matches the RingGrid flat index. function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractFullGrid) - polygons_matrix = RingGrids.get_gridcell_polygons(grid) - lonlat_to_usp = GO.UnitSphereFromGeographic() - + polygons = RingGrids.get_gridcell_polygons(grid) nlat = RingGrids.get_nlat(grid) nlon = RingGrids.get_nlon(grid) - poly_matrix = [_to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, (j - 1) * nlon + i) + poly_matrix = [ccw_unit_sphere_polygon(polygons, (j - 1) * nlon + i) for i in 1:nlon, j in 1:nlat] - epg = Trees.ExplicitPolygonGrid(manifold, poly_matrix) + epg = Trees.ExplicitPolygonGrid(manifold, poly_matrix) tree = Trees.TopDownQuadtreeCursor(epg) return Trees.KnownFullSphereExtentWrapper(tree) end @@ -42,14 +42,8 @@ end # wrapped in IndexOffsetQuadtreeCursor so that emitted indices match RingGrid flat indices. # ncells = npoints exactly — no ghost cells, no size mismatch in regrid!. function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) - polygons_matrix = RingGrids.get_gridcell_polygons(grid) - lonlat_to_usp = GO.UnitSphereFromGeographic() - - rings = RingGrids.eachring(grid) - nlat = length(rings) - - # Build one sub-tree per latitude ring - poly_type = typeof(_to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, 1)) + polygons = RingGrids.get_gridcell_polygons(grid) + rings = RingGrids.eachring(grid) subtrees = Trees.IndexOffsetQuadtreeCursor[] offsets = Int[] @@ -57,11 +51,9 @@ function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) for ring in rings nlon_ring = length(ring) - poly_ring = Matrix{poly_type}(undef, nlon_ring, 1) - for (i, ij) in enumerate(ring) - poly_ring[i, 1] = _to_ccw_unit_sphere_polygon(lonlat_to_usp, polygons_matrix, ij) - end - epg = Trees.ExplicitPolygonGrid(manifold, poly_ring) + poly_ring = [ccw_unit_sphere_polygon(polygons, ij) for ij in ring] + poly_ring = reshape(poly_ring, nlon_ring, 1) + epg = Trees.ExplicitPolygonGrid(manifold, poly_ring) cursor = Trees.IndexOffsetQuadtreeCursor(epg, first(ring) - 1) push!(subtrees, cursor) cumulative += nlon_ring From a2359a11014f024615bdd843b5c49e187e1fea52 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:27:23 +0200 Subject: [PATCH 07/13] Add GeoInterface, GeometryOps, GeometryOpsCore to deps Extensions can only use packages listed in the parent Project.toml. These are already in the dependency tree via ConservativeRegridding. Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 6 ++++++ .../speedy_regridder.jl | 14 +++++++------- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 7ce7d5615..3ad232040 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,9 @@ CFTime = "179af706-886a-5703-950a-314cd64e0468" ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" ConservativeRegridding = "8e50ac2c-eb48-49bc-a402-07c87b949343" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" +GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" +GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" +GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" @@ -81,6 +84,9 @@ StaticArrays = "1" Statistics = "<0.0.1, 1" Thermodynamics = "0.15.3" ConservativeRegridding = "0.2" +GeoInterface = "1" +GeometryOps = "0.1.33" +GeometryOpsCore = "0.1.9" WorldOceanAtlasTools = "0.6" ZipFile = "0.10" julia = "1.10" diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl index d9180b7a3..bd105f2de 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl @@ -6,8 +6,6 @@ import GeoInterface as GI using StaticArrays: SA -const lonlat_to_unit_sphere = GO.UnitSphereFromGeographic() - GOCore.best_manifold(grid::RingGrids.AbstractGrid) = GO.Spherical() GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.grid) @@ -15,10 +13,13 @@ treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(ma # get_gridcell_polygons returns CW (E, S, W, N); reverse to CCW (E, N, W, S) function ccw_unit_sphere_polygon(polygons_matrix, ij) - vE = lonlat_to_unit_sphere(polygons_matrix[1, ij]) - vN = lonlat_to_unit_sphere(polygons_matrix[4, ij]) - vW = lonlat_to_unit_sphere(polygons_matrix[3, ij]) - vS = lonlat_to_unit_sphere(polygons_matrix[2, ij]) + + USFG = GO.UnitSphereFromGeographic() + + vE = USFG(polygons_matrix[1, ij]) + vN = USFG(polygons_matrix[4, ij]) + vW = USFG(polygons_matrix[3, ij]) + vS = USFG(polygons_matrix[2, ij]) return GI.Polygon(SA[GI.LinearRing(SA[vE, vN, vW, vS, vE])]) end @@ -44,7 +45,6 @@ end function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) polygons = RingGrids.get_gridcell_polygons(grid) rings = RingGrids.eachring(grid) - subtrees = Trees.IndexOffsetQuadtreeCursor[] offsets = Int[] cumulative = 0 From 4d7e2b31355a447db175d991c855d5b81d1ac375 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 09:56:02 +0200 Subject: [PATCH 08/13] Add LinearAlgebra to deps for transpose in extension Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 3ad232040..b9e8a18bd 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" MeshArrays = "cb8c808f-1acf-59a3-9d2b-6e38d009f683" From 7cf0aae289e6e1476360b106a564a9e4936846b4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 14 Apr 2026 12:04:08 +0200 Subject: [PATCH 09/13] Fix radius mismatch between Oceananigans and SpeedyWeather manifolds Oceananigans uses radius=6.371e6, SpeedyWeather uses 6.3710088e6. Pass the exchange_grid's manifold explicitly to Regridder to avoid the "manifolds must be the same" error. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../speedy_weather_exchanger.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl index cc2b86ed7..8a700f695 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl @@ -27,7 +27,10 @@ net_fluxes(::SpeedySimulation) = nothing function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) spectral_grid = atmosphere.model.spectral_grid - from_atmosphere = Regridder(exchange_grid, spectral_grid) + # Use the exchange_grid's manifold for both grids to avoid + # radius mismatch between Oceananigans and SpeedyWeather. + manifold = GOCore.best_manifold(exchange_grid) + from_atmosphere = Regridder(manifold, exchange_grid, spectral_grid) to_atmosphere = transpose(from_atmosphere) regridder = (; to_atmosphere, from_atmosphere) From d25c51ec5b665a9e0bf9cc4044153b56dbde499e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 27 Apr 2026 10:36:09 +0200 Subject: [PATCH 10/13] Update spectral grid configuration in simulation --- examples/global_climate_simulation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl index b0fa5fbea..e8ae24f94 100644 --- a/examples/global_climate_simulation.jl +++ b/examples/global_climate_simulation.jl @@ -76,7 +76,7 @@ Oceananigans.set!(sea_ice.model, h=Metadatum(:sea_ice_thickness, dataset=ECCO4Mo # The `atmosphere_simulation` function takes care of building an atmosphere model with appropriate # hooks so that NumericalEarth can compute inter-component fluxes. nlayers = 4 -spectral_grid = SpeedyWeather.SpectralGrid(; trunc=63, nlayers) +spectral_grid = SpeedyWeather.SpectralGrid(; trunc=63, nlayers, Grid=FullClenshawGrid) atmosphere = atmosphere_simulation(spectral_grid, output=true) # The atmosphere model already includes some initial conditions as described above: From cba76b861533431884ea2da9c9409252fdeea331 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 1 May 2026 17:41:05 +0200 Subject: [PATCH 11/13] update the regridding --- Project.toml | 12 ++-- .../NumericalEarthSpeedyWeatherExt.jl | 1 - .../speedy_regridder.jl | 64 ------------------- 3 files changed, 4 insertions(+), 73 deletions(-) delete mode 100644 ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl diff --git a/Project.toml b/Project.toml index 8e55470bd..465741d10 100644 --- a/Project.toml +++ b/Project.toml @@ -11,9 +11,6 @@ ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" ConservativeRegridding = "8e50ac2c-eb48-49bc-a402-07c87b949343" CubedSphere = "7445602f-e544-4518-8976-18f8e8ae6cdb" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" -GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" -GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" -GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -21,7 +18,6 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" MeshArrays = "cb8c808f-1acf-59a3-9d2b-6e38d009f683" @@ -56,6 +52,9 @@ NumericalEarthSpeedyWeatherExt = "SpeedyWeather" NumericalEarthVerosExt = ["PythonCall", "CondaPkg"] NumericalEarthWOAExt = "WorldOceanAtlasTools" +[sources] +ConservativeRegridding = {url = "https://github.com/JuliaGeo/ConservativeRegridding.jl", rev = "main"} + [compat] Adapt = "4" Breeze = "0.4" @@ -88,10 +87,7 @@ SpeedyWeather = "0.19" StaticArrays = "1" Statistics = "<0.0.1, 1" Thermodynamics = "0.15.3" -ConservativeRegridding = "0.2" -GeoInterface = "1" -GeometryOps = "0.1.33" -GeometryOpsCore = "0.1.9" +ConservativeRegridding = "0.2.1" WorldOceanAtlasTools = "0.6" ZipFile = "0.10" julia = "1.10" diff --git a/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl b/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl index 047295a94..a980103a3 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/NumericalEarthSpeedyWeatherExt.jl @@ -10,7 +10,6 @@ import Oceananigans import SpeedyWeather.RingGrids include("speedy_atmosphere_simulations.jl") -include("speedy_regridder.jl") include("speedy_weather_exchanger.jl") end # module NumericalEarthSpeedyWeatherExt diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl deleted file mode 100644 index bd105f2de..000000000 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_regridder.jl +++ /dev/null @@ -1,64 +0,0 @@ -using ConservativeRegridding: Trees, Regridder -import ConservativeRegridding.Trees: treeify -import GeometryOpsCore as GOCore -import GeometryOps as GO -import GeoInterface as GI - -using StaticArrays: SA - -GOCore.best_manifold(grid::RingGrids.AbstractGrid) = GO.Spherical() -GOCore.best_manifold(sg::SpeedyWeather.SpectralGrid) = GOCore.best_manifold(sg.grid) - -treeify(manifold::GOCore.Spherical, sg::SpeedyWeather.SpectralGrid) = treeify(manifold, sg.grid) - -# get_gridcell_polygons returns CW (E, S, W, N); reverse to CCW (E, N, W, S) -function ccw_unit_sphere_polygon(polygons_matrix, ij) - - USFG = GO.UnitSphereFromGeographic() - - vE = USFG(polygons_matrix[1, ij]) - vN = USFG(polygons_matrix[4, ij]) - vW = USFG(polygons_matrix[3, ij]) - vS = USFG(polygons_matrix[2, ij]) - return GI.Polygon(SA[GI.LinearRing(SA[vE, vN, vW, vS, vE])]) -end - -# Full grids: all rings have the same nlon, so no padding needed. -# Linear index in ExplicitPolygonGrid matches the RingGrid flat index. -function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractFullGrid) - polygons = RingGrids.get_gridcell_polygons(grid) - nlat = RingGrids.get_nlat(grid) - nlon = RingGrids.get_nlon(grid) - - poly_matrix = [ccw_unit_sphere_polygon(polygons, (j - 1) * nlon + i) - for i in 1:nlon, j in 1:nlat] - - epg = Trees.ExplicitPolygonGrid(manifold, poly_matrix) - tree = Trees.TopDownQuadtreeCursor(epg) - return Trees.KnownFullSphereExtentWrapper(tree) -end - -# Reduced grids: variable nlon per ring. Build one sub-tree per latitude ring -# and combine with MultiTreeWrapper. Each ring is an ExplicitPolygonGrid(nlon_ring, 1) -# wrapped in IndexOffsetQuadtreeCursor so that emitted indices match RingGrid flat indices. -# ncells = npoints exactly — no ghost cells, no size mismatch in regrid!. -function treeify(manifold::GOCore.Spherical, grid::RingGrids.AbstractGrid) - polygons = RingGrids.get_gridcell_polygons(grid) - rings = RingGrids.eachring(grid) - subtrees = Trees.IndexOffsetQuadtreeCursor[] - offsets = Int[] - cumulative = 0 - - for ring in rings - nlon_ring = length(ring) - poly_ring = [ccw_unit_sphere_polygon(polygons, ij) for ij in ring] - poly_ring = reshape(poly_ring, nlon_ring, 1) - epg = Trees.ExplicitPolygonGrid(manifold, poly_ring) - cursor = Trees.IndexOffsetQuadtreeCursor(epg, first(ring) - 1) - push!(subtrees, cursor) - cumulative += nlon_ring - push!(offsets, cumulative) - end - - return Trees.MultiTreeWrapper(subtrees, offsets) -end From cf0c7b631e8af31014e75005f51b8b9ecb4bf3bd Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 1 May 2026 17:55:20 +0200 Subject: [PATCH 12/13] fix this --- .../speedy_weather_exchanger.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl index 3f4c8fb8d..9a7b698e6 100644 --- a/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl +++ b/ext/NumericalEarthSpeedyWeatherExt/speedy_weather_exchanger.jl @@ -5,8 +5,6 @@ using Oceananigans.Utils: launch! using Oceananigans.Operators: intrinsic_vector using ConservativeRegridding: Regridder -using LinearAlgebra: transpose - import ConservativeRegridding: regrid! using NumericalEarth.EarthSystemModels: sea_ice_concentration @@ -31,7 +29,7 @@ function ComponentExchanger(atmosphere::SpeedySimulation, exchange_grid) # radius mismatch between Oceananigans and SpeedyWeather. manifold = GOCore.best_manifold(exchange_grid) from_atmosphere = Regridder(manifold, exchange_grid, spectral_grid) - to_atmosphere = transpose(from_atmosphere) + to_atmosphere = Regridder(manifold, spectral_grid, exchange_grid) regridder = (; to_atmosphere, from_atmosphere) state = (; u = Field{Center, Center, Nothing}(exchange_grid), From b61b53d1b5ddc2bdd606f600e8428458c637d0eb Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 1 May 2026 17:56:35 +0200 Subject: [PATCH 13/13] Apply suggestion from @simone-silvestri --- test/test_speedy_coupling.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_speedy_coupling.jl b/test/test_speedy_coupling.jl index 7230449a2..dc4f8e652 100644 --- a/test/test_speedy_coupling.jl +++ b/test/test_speedy_coupling.jl @@ -7,7 +7,7 @@ using Test NumericalEarthSpeedyWeatherExt = Base.get_extension(NumericalEarth, :NumericalEarthSpeedyWeatherExt) @test !isnothing(NumericalEarthSpeedyWeatherExt) -spectral_grid = SpeedyWeather.SpectralGrid(trunc=51, nlayers=3) +spectral_grid = SpeedyWeather.SpectralGrid(trunc=51, nlayers=3, Grid=FullClenshawGrid) oceananigans_grid = LatitudeLongitudeGrid(Oceananigans.CPU(); size=(200, 100, 1), latitude=(-80, 80), longitude=(0, 360), z = (0, 1)) ocean = NumericalEarth.Oceans.ocean_simulation(oceananigans_grid; momentum_advection=nothing, tracer_advection=nothing, closure=nothing)