From 3e1178b21fcb74ca6b24937dc0d0c75cb19bdcd1 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 14:14:47 +0100 Subject: [PATCH 01/22] one missing allowscalar --- .../InterfaceComputations/component_interfaces.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index f83e0c84..3bbe5cc0 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -83,7 +83,7 @@ function StateExchanger(ocean::Simulation, atmosphere) # Make an array of FractionalIndices kᴺ = size(exchange_grid, 3) @allowscalar X1 = _node(1, 1, kᴺ + 1, exchange_grid, c, c, f) - i1 = FractionalIndices(X1, atmosphere.grid, c, c, nothing) + @allowscalar i1 = FractionalIndices(X1, atmosphere.grid, c, c, nothing) frac_indices_data = [deepcopy(i1) for i=1:Nx+2, j=1:Ny+2, k=1:1] frac_indices = OffsetArray(frac_indices_data, -1, -1, 0) frac_indices = on_architecture(arch, frac_indices) From 72dcc2fc34f7043cd43c049ff8359cb1df2cc884 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 14:25:46 +0100 Subject: [PATCH 02/22] small cleanup --- .../InterfaceComputations/component_interfaces.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 3bbe5cc0..91057964 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -2,7 +2,6 @@ using StaticArrays using Thermodynamics using SurfaceFluxes using OffsetArrays -using CUDA: @allowscalar using ..OceanSeaIceModels: reference_density, heat_capacity, @@ -81,10 +80,8 @@ function StateExchanger(ocean::Simulation, atmosphere) Nx, Ny, Nz = size(exchange_grid) # Make an array of FractionalIndices - kᴺ = size(exchange_grid, 3) - @allowscalar X1 = _node(1, 1, kᴺ + 1, exchange_grid, c, c, f) - @allowscalar i1 = FractionalIndices(X1, atmosphere.grid, c, c, nothing) - frac_indices_data = [deepcopy(i1) for i=1:Nx+2, j=1:Ny+2, k=1:1] + FT = eltype(exchange_grid) + frac_indices = [FractionalIndices(one(FT), one(FT), one(FT)) for i=1:Nx+2, j=1:Ny+2, k=1:1] frac_indices = OffsetArray(frac_indices_data, -1, -1, 0) frac_indices = on_architecture(arch, frac_indices) From b0ef88fe83ce0de0e08c7c1f5f051abc49603415 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 14:28:48 +0100 Subject: [PATCH 03/22] bugfix --- .../InterfaceComputations/component_interfaces.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 91057964..689e17b9 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -82,7 +82,7 @@ function StateExchanger(ocean::Simulation, atmosphere) # Make an array of FractionalIndices FT = eltype(exchange_grid) frac_indices = [FractionalIndices(one(FT), one(FT), one(FT)) for i=1:Nx+2, j=1:Ny+2, k=1:1] - frac_indices = OffsetArray(frac_indices_data, -1, -1, 0) + frac_indices = OffsetArray(frac_indices, -1, -1, 0) frac_indices = on_architecture(arch, frac_indices) kernel_parameters = interface_kernel_parameters(exchange_grid) From 3d07bc70200db917652f5b5a0a118de265db5ba0 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 4 Mar 2025 07:19:58 -0700 Subject: [PATCH 04/22] New way of constructing fractional indices --- .../InterfaceComputations/component_interfaces.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 689e17b9..20b70276 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -59,6 +59,11 @@ struct StateExchanger{G, AST, AEX} atmosphere_exchanger :: AEX end +# Note that Field location can also affect fractional index type. +# Here we assume that we know the location of Fields that will be interpolated. +fractional_index_type(FT, Topo) = FT +fractional_index_type(FT, ::Flat) = Nothing + function StateExchanger(ocean::Simulation, atmosphere) # TODO: generalize this exchange_grid = ocean.model.grid @@ -81,8 +86,13 @@ function StateExchanger(ocean::Simulation, atmosphere) # Make an array of FractionalIndices FT = eltype(exchange_grid) - frac_indices = [FractionalIndices(one(FT), one(FT), one(FT)) for i=1:Nx+2, j=1:Ny+2, k=1:1] - frac_indices = OffsetArray(frac_indices, -1, -1, 0) + TX, TY, TZ = topology(exchange_grid) + FX = fractional_index_type(FT, TX()) + FY = fractional_index_type(FT, TY()) + FR = FractionalIndices{FX, FY, Nothing} + + underlying_frac_indices = Array{FR}(undef, Nx+2, Ny+2, 1) + frac_indices = OffsetArray(underlying_frac_indices, -1, -1, 0) frac_indices = on_architecture(arch, frac_indices) kernel_parameters = interface_kernel_parameters(exchange_grid) From cf656fe4a186dcfbf4f62a370fb0f7de2e54bc17 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 16:12:45 +0100 Subject: [PATCH 05/22] add topoology --- .../InterfaceComputations/component_interfaces.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 20b70276..ac30239d 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -14,7 +14,7 @@ using ..OceanSeaIceModels: reference_density, using ClimaSeaIce: SeaIceModel using Oceananigans: HydrostaticFreeSurfaceModel, architecture -using Oceananigans.Grids: inactive_node, node +using Oceananigans.Grids: inactive_node, node, topology using Oceananigans.BoundaryConditions: fill_halo_regions! using Oceananigans.Fields: ConstantField, interpolate, FractionalIndices From c4f1e932e4cd923ea80e779f0f3a8dc7f3167c5b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 4 Mar 2025 08:30:34 -0700 Subject: [PATCH 06/22] Tuple of fractional indices --- .../component_interfaces.jl | 24 +++++++----- .../interpolate_atmospheric_state.jl | 38 +++++++++++-------- 2 files changed, 37 insertions(+), 25 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index ac30239d..f9096ddf 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -84,16 +84,16 @@ function StateExchanger(ocean::Simulation, atmosphere) arch = architecture(exchange_grid) Nx, Ny, Nz = size(exchange_grid) - # Make an array of FractionalIndices - FT = eltype(exchange_grid) + # Make a NamedTuple of fractional indices + # Note: we could use an array of FractionalIndices. Instead, for compatbility + # with Reactant we construct FractionalIndices on the fly in `interpolate_atmospheric_state`. + FT = eltype(atmos_grid) TX, TY, TZ = topology(exchange_grid) FX = fractional_index_type(FT, TX()) FY = fractional_index_type(FT, TY()) - FR = FractionalIndices{FX, FY, Nothing} - - underlying_frac_indices = Array{FR}(undef, Nx+2, Ny+2, 1) - frac_indices = OffsetArray(underlying_frac_indices, -1, -1, 0) - frac_indices = on_architecture(arch, frac_indices) + fi = Field{Center, Center, Nothing}(exchange_grid, FX) + fj = Field{Center, Center, Nothing}(exchange_grid, FY) + frac_indices = (i=fi, j=fj) # no k needed, only horizontal interpolation kernel_parameters = interface_kernel_parameters(exchange_grid) launch!(arch, exchange_grid, kernel_parameters, @@ -102,11 +102,17 @@ function StateExchanger(ocean::Simulation, atmosphere) return StateExchanger(ocean.model.grid, exchange_atmosphere_state, frac_indices) end -@kernel function _compute_fractional_indices!(frac_indices, exchange_grid, atmos_grid) +@kernel function _compute_fractional_indices!(indices_tuple, exchange_grid, atmos_grid) i, j = @index(Global, NTuple) kᴺ = size(exchange_grid, 3) # index of the top ocean cell X = _node(i, j, kᴺ + 1, exchange_grid, c, c, f) - @inbounds frac_indices[i, j, 1] = FractionalIndices(X, atmos_grid, c, c, nothing) + fractional_indices_ij = FractionalIndices(X, atmos_grid, c, c, nothing) + fi = indices_tuple.i + fj = indices_tuple.j + @inbounds begin + fi[i, j, k] = fractional_indices_ij.i + fj[i, j, k] = fractional_indices_ij.j + end end const PATP = PrescribedAtmosphereThermodynamicsParameters diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index 3bc69250..9c800054 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -1,5 +1,6 @@ using Oceananigans.Operators: intrinsic_vector using Oceananigans.Grids: _node +using Oceananigans.Fields: FractionalIndices using Oceananigans.OutputReaders: TimeInterpolator using ...OceanSimulations: forcing_barotropic_potential @@ -125,27 +126,32 @@ end i, j = @index(Global, NTuple) @inbounds begin - x_itp = space_fractional_indices[i, j, 1] - t_itp = time_interpolator - atmos_args = (x_itp, t_itp, atmos_backend, atmos_time_indexing) + fi = space_fractional_indices.fi[i, j, 1] + fj = space_fractional_indices.fj[i, j, 1] + end + + x_itp = FractionalIndices(fi, fj, nothing) + t_itp = time_interpolator + atmos_args = (x_itp, t_itp, atmos_backend, atmos_time_indexing) - uₐ = interp_atmos_time_series(atmos_velocities.u, atmos_args...) - vₐ = interp_atmos_time_series(atmos_velocities.v, atmos_args...) - Tₐ = interp_atmos_time_series(atmos_tracers.T, atmos_args...) - qₐ = interp_atmos_time_series(atmos_tracers.q, atmos_args...) - pₐ = interp_atmos_time_series(atmos_pressure, atmos_args...) + uₐ = interp_atmos_time_series(atmos_velocities.u, atmos_args...) + vₐ = interp_atmos_time_series(atmos_velocities.v, atmos_args...) + Tₐ = interp_atmos_time_series(atmos_tracers.T, atmos_args...) + qₐ = interp_atmos_time_series(atmos_tracers.q, atmos_args...) + pₐ = interp_atmos_time_series(atmos_pressure, atmos_args...) - Qs = interp_atmos_time_series(downwelling_radiation.shortwave, atmos_args...) - Qℓ = interp_atmos_time_series(downwelling_radiation.longwave, atmos_args...) + Qs = interp_atmos_time_series(downwelling_radiation.shortwave, atmos_args...) + Qℓ = interp_atmos_time_series(downwelling_radiation.longwave, atmos_args...) - # Usually precipitation - Mh = interp_atmos_time_series(prescribed_freshwater_flux, atmos_args...) + # Usually precipitation + Mh = interp_atmos_time_series(prescribed_freshwater_flux, atmos_args...) - # Convert atmosphere velocities (usually defined on a latitude-longitude grid) to - # the frame of reference of the native grid - kᴺ = size(exchange_grid, 3) # index of the top ocean cell - uₐ, vₐ = intrinsic_vector(i, j, kᴺ, exchange_grid, uₐ, vₐ) + # Convert atmosphere velocities (usually defined on a latitude-longitude grid) to + # the frame of reference of the native grid + kᴺ = size(exchange_grid, 3) # index of the top ocean cell + uₐ, vₐ = intrinsic_vector(i, j, kᴺ, exchange_grid, uₐ, vₐ) + @inbounds begin surface_atmos_state.u[i, j, 1] = uₐ surface_atmos_state.v[i, j, 1] = vₐ surface_atmos_state.T[i, j, 1] = Tₐ From db2b3dc84bb98f0d9057ca4d94d4002da70924c2 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 4 Mar 2025 12:15:18 -0700 Subject: [PATCH 07/22] Bugfix --- .../InterfaceComputations/component_interfaces.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index f9096ddf..26347cf6 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -110,8 +110,8 @@ end fi = indices_tuple.i fj = indices_tuple.j @inbounds begin - fi[i, j, k] = fractional_indices_ij.i - fj[i, j, k] = fractional_indices_ij.j + fi[i, j, 1] = fractional_indices_ij.i + fj[i, j, 1] = fractional_indices_ij.j end end From a2f35f7938fdd6a5222e3dee599065d32f79bd00 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 4 Mar 2025 14:52:31 -0500 Subject: [PATCH 08/22] Bump ClimaSeaICe --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 134ee327..66b101e6 100644 --- a/Project.toml +++ b/Project.toml @@ -33,7 +33,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" Adapt = "4" CFTime = "0.1" CUDA = "4, 5" -ClimaSeaIce = "0.2.0" +ClimaSeaIce = "0.2.1" CubicSplines = "0.2" DataDeps = "0.7" Downloads = "1.6" From 9016c17335b27bebb362d6cbe4576b890fe8cd1a Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 4 Mar 2025 14:52:44 -0500 Subject: [PATCH 09/22] Fix bug in atmos interpolation --- .../InterfaceComputations/interpolate_atmospheric_state.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index 9c800054..5d7366ba 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -126,8 +126,8 @@ end i, j = @index(Global, NTuple) @inbounds begin - fi = space_fractional_indices.fi[i, j, 1] - fj = space_fractional_indices.fj[i, j, 1] + fi = space_fractional_indices.i[i, j, 1] + fj = space_fractional_indices.j[i, j, 1] end x_itp = FractionalIndices(fi, fj, nothing) From 788aa896edc3bd7f6681adc2b22c3361bb30b85b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 4 Mar 2025 19:00:15 -0500 Subject: [PATCH 10/22] Try to generalize to doubly-Flat grids --- .../InterfaceComputations/component_interfaces.jl | 15 +++++++++------ .../interpolate_atmospheric_state.jl | 11 +++++++---- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 26347cf6..0f74f6ee 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -89,10 +89,8 @@ function StateExchanger(ocean::Simulation, atmosphere) # with Reactant we construct FractionalIndices on the fly in `interpolate_atmospheric_state`. FT = eltype(atmos_grid) TX, TY, TZ = topology(exchange_grid) - FX = fractional_index_type(FT, TX()) - FY = fractional_index_type(FT, TY()) - fi = Field{Center, Center, Nothing}(exchange_grid, FX) - fj = Field{Center, Center, Nothing}(exchange_grid, FY) + fi = TX() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FX) + fj = TY() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FY) frac_indices = (i=fi, j=fj) # no k needed, only horizontal interpolation kernel_parameters = interface_kernel_parameters(exchange_grid) @@ -110,8 +108,13 @@ end fi = indices_tuple.i fj = indices_tuple.j @inbounds begin - fi[i, j, 1] = fractional_indices_ij.i - fj[i, j, 1] = fractional_indices_ij.j + if !isnothing(fi) + fi[i, j, 1] = fractional_indices_ij.i + end + + if !isnothing(fj) + fj[i, j, 1] = fractional_indices_ij.j + end end end diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index 5d7366ba..193276a9 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -110,6 +110,9 @@ function interpolate_atmospheric_state!(coupled_model) parent(barotropic_potential) .= parent(atmosphere_data.p) ./ ρₒ end end + +@inline get_fractional_index(i, j, ::Nothing) = 1 +@inline get_fractional_index(i, j, frac) = @inbounds frac[i, j, 1] @kernel function _interpolate_primary_atmospheric_state!(surface_atmos_state, space_fractional_indices, @@ -125,10 +128,10 @@ end i, j = @index(Global, NTuple) - @inbounds begin - fi = space_fractional_indices.i[i, j, 1] - fj = space_fractional_indices.j[i, j, 1] - end + ii = space_fractional_indices.i + jj = space_fractional_indices.j + fi = get_fractional_index(i, j, ii) + fj = get_fractional_index(i, j, jj) x_itp = FractionalIndices(fi, fj, nothing) t_itp = time_interpolator From 952e901e61db98c274fe25295f259e1406c569f9 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 09:54:00 -0700 Subject: [PATCH 11/22] Bug in StateExchanger construcotr --- .../InterfaceComputations/component_interfaces.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 0f74f6ee..097e6bb1 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -89,8 +89,8 @@ function StateExchanger(ocean::Simulation, atmosphere) # with Reactant we construct FractionalIndices on the fly in `interpolate_atmospheric_state`. FT = eltype(atmos_grid) TX, TY, TZ = topology(exchange_grid) - fi = TX() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FX) - fj = TY() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FY) + fi = TX() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FT) + fj = TY() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FT) frac_indices = (i=fi, j=fj) # no k needed, only horizontal interpolation kernel_parameters = interface_kernel_parameters(exchange_grid) From ffd937309358b525c75e359b381694ce16f456d3 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 09:54:23 -0700 Subject: [PATCH 12/22] Bump Oceananigans compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 66b101e6..006bdad4 100644 --- a/Project.toml +++ b/Project.toml @@ -42,7 +42,7 @@ JLD2 = "0.4, 0.5" KernelAbstractions = "0.9" MPI = "0.20" NCDatasets = "0.12, 0.13, 0.14" -Oceananigans = "0.95.17 - 0.99" +Oceananigans = "0.95.18 - 0.99" OffsetArrays = "1.14" OrthogonalSphericalShellGrids = "0.2.2" Scratch = "1" From 360bfa6dfe30d1ba441b44a87477259830b8f81b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 5 Mar 2025 20:37:37 +0100 Subject: [PATCH 13/22] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 006bdad4..68dc0fd3 100644 --- a/Project.toml +++ b/Project.toml @@ -33,7 +33,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" Adapt = "4" CFTime = "0.1" CUDA = "4, 5" -ClimaSeaIce = "0.2.1" +ClimaSeaIce = "0.2.0" CubicSplines = "0.2" DataDeps = "0.7" Downloads = "1.6" From 635ead9494b9a48fb9073246f7ecb7e40025f8b2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 5 Mar 2025 21:51:45 +0100 Subject: [PATCH 14/22] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 68dc0fd3..006bdad4 100644 --- a/Project.toml +++ b/Project.toml @@ -33,7 +33,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" Adapt = "4" CFTime = "0.1" CUDA = "4, 5" -ClimaSeaIce = "0.2.0" +ClimaSeaIce = "0.2.1" CubicSplines = "0.2" DataDeps = "0.7" Downloads = "1.6" From e4822b577fd3734077c3f0838a8728c6cbc7f8d6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 5 Mar 2025 21:53:42 +0100 Subject: [PATCH 15/22] fix test simulation --- test/test_simulations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index c0b2f599..8dbe5e06 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -56,8 +56,8 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature above_freezing_ocean_temperature!(ocean, sea_ice) # Test that ocean temperatures are above freezing - T = on_architecture(CPU(), ocean.model.T) - S = on_architecture(CPU(), ocean.model.S) + T = on_architecture(CPU(), ocean.model.tracers.T) + S = on_architecture(CPU(), ocean.model.tracers.S) @inline pointwise_melting_T(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) From aa0836026736bf56fbd7d06060369a372c84b4f5 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 14:36:28 -0700 Subject: [PATCH 16/22] Fix interpolation for Flat topologies? --- .../InterfaceComputations/interpolate_atmospheric_state.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index 193276a9..c027af12 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -111,7 +111,7 @@ function interpolate_atmospheric_state!(coupled_model) end end -@inline get_fractional_index(i, j, ::Nothing) = 1 +@inline get_fractional_index(i, j, ::Nothing) = 0 @inline get_fractional_index(i, j, frac) = @inbounds frac[i, j, 1] @kernel function _interpolate_primary_atmospheric_state!(surface_atmos_state, From b6e57c4e63cc1650c7f2f5fe3fe0639cb7f8aa68 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 15:28:24 -0700 Subject: [PATCH 17/22] Ohhh --- .../InterfaceComputations/interpolate_atmospheric_state.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index c027af12..308ad691 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -111,7 +111,7 @@ function interpolate_atmospheric_state!(coupled_model) end end -@inline get_fractional_index(i, j, ::Nothing) = 0 +@inline get_fractional_index(i, j, ::Nothing) = nothing @inline get_fractional_index(i, j, frac) = @inbounds frac[i, j, 1] @kernel function _interpolate_primary_atmospheric_state!(surface_atmos_state, From b3b3b30255dd2af8254582b2e87c0ed98d2c7f96 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 16:36:52 -0700 Subject: [PATCH 18/22] Fix test bug --- test/test_simulations.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 8dbe5e06..5b5d0f81 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -5,6 +5,8 @@ using OrthogonalSphericalShellGrids using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature! using ClimaSeaIce.SeaIceThermodynamics: melting_temperature +@inline pointwise_melting_T(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) + @testset "GPU time stepping test" begin for arch in test_architectures @@ -59,10 +61,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature T = on_architecture(CPU(), ocean.model.tracers.T) S = on_architecture(CPU(), ocean.model.tracers.S) - @inline pointwise_melting_T(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) - - Tm = KernelFunctionOperation{Center, Center, Center}(pointwise_melting_T, grid, S) - + Tm = KernelFunctionOperation{Center, Center, Center}(pointwise_melting_T, grid, liquidus, S) @test all(T .> Tm) # Fluxes are computed when the model is constructed, so we just test that this works. From 2b57bb1de96d3bf97ece78e1781daccd1b205fe1 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 16:37:35 -0700 Subject: [PATCH 19/22] Better name --- test/test_simulations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 5b5d0f81..e12c7771 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -5,7 +5,7 @@ using OrthogonalSphericalShellGrids using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature! using ClimaSeaIce.SeaIceThermodynamics: melting_temperature -@inline pointwise_melting_T(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) +@inline melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) @testset "GPU time stepping test" begin @@ -61,7 +61,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature T = on_architecture(CPU(), ocean.model.tracers.T) S = on_architecture(CPU(), ocean.model.tracers.S) - Tm = KernelFunctionOperation{Center, Center, Center}(pointwise_melting_T, grid, liquidus, S) + Tm = KernelFunctionOperation{Center, Center, Center}(melting_temperature, grid, liquidus, S) @test all(T .> Tm) # Fluxes are computed when the model is constructed, so we just test that this works. From 85f641e7972f43a21754749c02a36c3c08a68a8c Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 16:38:53 -0700 Subject: [PATCH 20/22] Name mangle --- test/test_simulations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index e12c7771..e7c41720 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -5,7 +5,7 @@ using OrthogonalSphericalShellGrids using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature! using ClimaSeaIce.SeaIceThermodynamics: melting_temperature -@inline melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) +@inline kernel_melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) @testset "GPU time stepping test" begin @@ -61,7 +61,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature T = on_architecture(CPU(), ocean.model.tracers.T) S = on_architecture(CPU(), ocean.model.tracers.S) - Tm = KernelFunctionOperation{Center, Center, Center}(melting_temperature, grid, liquidus, S) + Tm = KernelFunctionOperation{Center, Center, Center}(kernel_melting_temperature, grid, liquidus, S) @test all(T .> Tm) # Fluxes are computed when the model is constructed, so we just test that this works. From 990deb95708f1a39d04aebea66b85634316fe5b2 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 16:53:19 -0700 Subject: [PATCH 21/22] Fix test --- test/test_simulations.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index e7c41720..75681f51 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -7,7 +7,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature @inline kernel_melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) -@testset "GPU time stepping test" begin +@testset "Time stepping test" begin for arch in test_architectures @@ -26,9 +26,9 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature interpolation_passes = 20, major_basins = 1) - grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true) + grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) - free_surface = SplitExplicitFreeSurface(grid; substeps = 20) + free_surface = SplitExplicitFreeSurface(grid; substeps=20) ocean = ocean_simulation(grid; free_surface) backend = JRA55NetCDFBackend(4) @@ -54,7 +54,6 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature # Set the ocean temperature and salinity set!(ocean.model, T=temperature_metadata[1], S=salinity_metadata[1]) - above_freezing_ocean_temperature!(ocean, sea_ice) # Test that ocean temperatures are above freezing @@ -62,13 +61,15 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature S = on_architecture(CPU(), ocean.model.tracers.S) Tm = KernelFunctionOperation{Center, Center, Center}(kernel_melting_temperature, grid, liquidus, S) - @test all(T .> Tm) + @test all(T .>= Tm) # Fluxes are computed when the model is constructed, so we just test that this works. # And that we can time step with sea ice @test begin coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + time_step!(coupled_model, 1) true end end end + From 9c0cc8998efbb554af748d6cf707d48de3f0258b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 5 Mar 2025 22:20:44 -0700 Subject: [PATCH 22/22] Tiny change in MixedLayerDepth diagnostic --- test/test_diagnostics.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_diagnostics.jl b/test/test_diagnostics.jl index fc64742d..ec9ef89a 100644 --- a/test/test_diagnostics.jl +++ b/test/test_diagnostics.jl @@ -42,11 +42,11 @@ using ClimaOcean.Diagnostics: MixedLayerDepthField, MixedLayerDepthOperand @test h.operand.buoyancy_perturbation isa KernelFunctionOperation compute!(h) - @test @allowscalar h[1, 1, 1] ≈ 16.255836 # m + @test @allowscalar h[1, 1, 1] ≈ 16.2558363 # m tracers = (T=Tt[2], S=St[2]) h.operand.buoyancy_perturbation = buoyancy(sb, grid, tracers) compute!(h) - @test @allowscalar h[1, 1, 1] ≈ 9.295890287 # m + @test @allowscalar h[1, 1, 1] ≈ 9.2957298 # m end end