From 922a36e900cc3724ecadead12fd7c1562a9df8b8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 31 Jan 2025 09:34:33 +0100 Subject: [PATCH 01/56] Update Bathymetry.jl --- src/Bathymetry.jl | 24 ++++-------------------- 1 file changed, 4 insertions(+), 20 deletions(-) diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index b9148d06a..42f884a1a 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -125,25 +125,9 @@ function regrid_bathymetry(target_grid; φ₁, φ₂ = y_domain(target_grid) λ₁, λ₂ = x_domain(target_grid) - if λ₁ < 0 || λ₂ > 360 - throw(ArgumentError("Cannot regrid bathymetry between λ₁ = $(λ₁) and λ₂ = $(λ₂). - Bathymetry data is defined on longitudes spanning λ = (0, 360).")) - end - - # Calculate limiting indices on the bathymetry grid - i₁ = searchsortedfirst(λ_data, λ₁) - i₂ = searchsortedfirst(λ_data, λ₂) - 1 - ii = i₁:i₂ - - j₁ = searchsortedfirst(φ_data, φ₁) - j₂ = searchsortedfirst(φ_data, φ₂) - 1 - jj = j₁:j₂ - - # Restrict bathymetry coordinate_data to region of interest - λ_data = λ_data[ii] |> Array{BigFloat} - φ_data = φ_data[jj] |> Array{BigFloat} - z_data = z_data[ii, jj] - + λ_data = λ_data |> Array{BigFloat} + φ_data = φ_data |> Array{BigFloat} + if !isnothing(height_above_water) # Overwrite the height of cells above water. # This has an impact on reconstruction. Greater height_above_water reduces total @@ -164,7 +148,7 @@ function regrid_bathymetry(target_grid; Nxn = length(λ_data) Nyn = length(φ_data) Nzn = 1 - + native_grid = LatitudeLongitudeGrid(arch; size = (Nxn, Nyn, Nzn), latitude = (φ₁_data, φ₂_data), From cc056e0b02a4d864df6aae1dd6803f209e073ad5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Feb 2025 10:30:23 +0100 Subject: [PATCH 02/56] add a fill halo --- src/Bathymetry.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index 42f884a1a..9b0ecd397 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -158,7 +158,8 @@ function regrid_bathymetry(target_grid; native_z = Field{Center, Center, Nothing}(native_grid) set!(native_z, z_data) - + fill_halo_regions!(native_z) + target_z = interpolate_bathymetry_in_passes(native_z, target_grid; passes = interpolation_passes) From e17590640cff046be9274fe5060185d20dc46dea Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 11:26:34 +0100 Subject: [PATCH 03/56] starting --- experiments/arctic_simulation.jl | 52 ++++++++++++++++++++++++++++++++ src/Bathymetry.jl | 6 ++-- 2 files changed, 54 insertions(+), 4 deletions(-) create mode 100644 experiments/arctic_simulation.jl diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl new file mode 100644 index 000000000..30541816d --- /dev/null +++ b/experiments/arctic_simulation.jl @@ -0,0 +1,52 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.OrthogonalSphericalShellGrids + +r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=3000) +z_faces = MutableVerticalDiscretization(r_faces) + +Nx = 180 # longitudinal direction -> 250 points is about 1.5ᵒ resolution +Ny = 180 # meridional direction -> same thing, 48 points is about 1.5ᵒ resolution +Nz = length(r_faces) - 1 + +grid = RotatedLatitudeLongitudeGrid(size = (Nx, Ny, Nz), + latitude = (-45, 45), + longitude = (-45, 45), + z = z_faces, + north_pole = (0, 0), + topology = (Bounded, Bounded, Bounded)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) + +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) + +@inline x_domain(grid::RotatedLatitudeLongitudeGrid) = domain(topology(grid, 1)(), grid.Nx, grid.xᶠᵃᵃ) +@inline y_domain(grid::RotatedLatitudeLongitudeGrid) = domain(topology(grid, 2)(), grid.Ny, grid.yᵃᶠᵃ) + +##### +##### Ocean model +##### + +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = Centered() + +free_surface = SplitExplicitFreeSurface(grid; substeps=30) + +ocean = ocean_simulation(grid; + momentum_advection, + tracer_advection, + free_surface) + +##### +##### Sea-ice model +##### + +sea_ice = sea_ice_simulation(grid; dynamics=nothing, advection=nothing) + +##### +##### Atmosphere model +##### + +atmosphere = JRA55PrescribedAtmosphere(; backend=JRA55NetCDFBackend(40)) diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index d109b6306..3ddaf5abd 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -127,8 +127,6 @@ function regrid_bathymetry(target_grid; # Diagnose target grid information arch = architecture(target_grid) - φ₁, φ₂ = y_domain(target_grid) - λ₁, λ₂ = x_domain(target_grid) λ_data = λ_data |> Array{BigFloat} φ_data = φ_data |> Array{BigFloat} @@ -205,8 +203,8 @@ function interpolate_bathymetry_in_passes(native_z, target_grid; end # Interpolate in passes - latitude = y_domain(target_grid) - longitude = x_domain(target_grid) + latitude = y_domain(native_z.grid) + longitude = x_domain(native_z.grid) ΔNλ = floor((Nλn - Nλt) / passes) ΔNφ = floor((Nφn - Nφt) / passes) From 96aed5b1a68723df3d34124a80126fb26bc16639 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 11:55:43 +0100 Subject: [PATCH 04/56] some changes --- experiments/arctic_simulation.jl | 28 +++++++++++++++++++++++----- src/SeaIceSimulations.jl | 1 + 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 30541816d..fc0b9d873 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -2,6 +2,7 @@ using ClimaOcean using ClimaSeaIce using Oceananigans using Oceananigans.Grids +using Oceananigans.Units using Oceananigans.OrthogonalSphericalShellGrids r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=3000) @@ -15,16 +16,14 @@ grid = RotatedLatitudeLongitudeGrid(size = (Nx, Ny, Nz), latitude = (-45, 45), longitude = (-45, 45), z = z_faces, - north_pole = (0, 0), + north_pole = (180, 0), + halo = (5, 5, 4), topology = (Bounded, Bounded, Bounded)) bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) -@inline x_domain(grid::RotatedLatitudeLongitudeGrid) = domain(topology(grid, 1)(), grid.Nx, grid.xᶠᵃᵃ) -@inline y_domain(grid::RotatedLatitudeLongitudeGrid) = domain(topology(grid, 2)(), grid.Ny, grid.yᵃᶠᵃ) - ##### ##### Ocean model ##### @@ -32,7 +31,7 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) momentum_advection = WENOVectorInvariant(order=5) tracer_advection = Centered() -free_surface = SplitExplicitFreeSurface(grid; substeps=30) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) ocean = ocean_simulation(grid; momentum_advection, @@ -45,8 +44,27 @@ ocean = ocean_simulation(grid; sea_ice = sea_ice_simulation(grid; dynamics=nothing, advection=nothing) +##### +##### Initialize the models +##### + +set!(ocean.model, T=ECCOMetadata(:temperature), + S=ECCOMetadata(:salinity)) + +set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), + ℵ=ECCOMetadata(:sea_ice_concentration)) + ##### ##### Atmosphere model ##### atmosphere = JRA55PrescribedAtmosphere(; backend=JRA55NetCDFBackend(40)) +radiation = Radiation() + +##### +##### Arctic coupled model +##### + +arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +arctic = Simulation(arctic, Δt=600, stop_time=365days) + diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index c46f5d3dc..834be7db1 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -16,6 +16,7 @@ using Oceananigans.Operators using ClimaSeaIce using ClimaSeaIce: SeaIceModel, SlabSeaIceThermodynamics, PhaseTransitions, ConductiveFlux +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium using ClimaOcean.OceanSimulations: Default From c89d8f9c47ea4623cb35c070e248094f5938ff46 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 12:24:32 +0100 Subject: [PATCH 05/56] start changing stuff --- experiments/arctic_simulation.jl | 73 +++++++++++++++---- .../flux_climatology/flux_climatology.jl | 57 +-------------- src/Bathymetry.jl | 4 +- src/OceanSimulations/OceanSimulations.jl | 2 +- src/OceanSimulations/prescribed_ocean.jl | 71 ++++++++++++++++++ 5 files changed, 133 insertions(+), 74 deletions(-) create mode 100644 src/OceanSimulations/prescribed_ocean.jl diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index fc0b9d873..d291e38b5 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -4,6 +4,10 @@ using Oceananigans using Oceananigans.Grids using Oceananigans.Units using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.ECCO: all_ECCO_dates +using Printf r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=3000) z_faces = MutableVerticalDiscretization(r_faces) @@ -25,32 +29,27 @@ bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) ##### -##### Ocean model +##### A Prescribed Ocean model ##### -momentum_advection = WENOVectorInvariant(order=5) -tracer_advection = Centered() +# ...with prescribed velocity and tracer fields +version = ECCO2Daily() +dates = all_ECCO_dates(version)[1:200] -free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) +u = ECCOFieldTimeSeries(:u_velocity, version; dates, time_indices_in_memory=2) +v = ECCOFieldTimeSeries(:v_velocity, version; dates, time_indices_in_memory=2) +T = ECCOFieldTimeSeries(:temperature, version; dates, time_indices_in_memory=2) +S = ECCOFieldTimeSeries(:salinity, version; dates, time_indices_in_memory=2) -ocean = ocean_simulation(grid; - momentum_advection, - tracer_advection, - free_surface) +ocean_model = PrescribedOcean((; u, v, T, S); grid) +ocean = Simulation(ocean_model, Δt=10minutes, verbose=false) ##### -##### Sea-ice model +##### A Prognostic Sea-ice model ##### sea_ice = sea_ice_simulation(grid; dynamics=nothing, advection=nothing) -##### -##### Initialize the models -##### - -set!(ocean.model, T=ECCOMetadata(:temperature), - S=ECCOMetadata(:salinity)) - set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), ℵ=ECCOMetadata(:sea_ice_concentration)) @@ -68,3 +67,45 @@ radiation = Radiation() arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) arctic = Simulation(arctic, Δt=600, stop_time=365days) +h = sea_ice.model.ice_thickness +ℵ = sea_ice.model.ice_concentration +Gh = sea_ice.model.timestepper.Gⁿ.h +Gℵ = sea_ice.model.timestepper.Gⁿ.ℵ +Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature + +sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu), + filename = "sea_ice_quantities.jld2", + schedule = IterationInterval(12)) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + hmean = mean(sea_ice.model.ice_thickness) + ℵmean = mean(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg3 = @sprintf("mean(h): %.2e m, mean(ℵ): %.2e ", hmean, ℵmean) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 * msg4 * msg5 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(arctic, progress, IterationInterval(10)) + +run!(arctic) \ No newline at end of file diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 94a1a951e..aff118f56 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -13,7 +13,7 @@ using Adapt using Printf using KernelAbstractions: @index, @kernel -import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! +import Oceananigans.TimeSteppers: time_step!, update_state! import Oceananigans.Models: timestepper, update_model_field_time_series! import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity @@ -115,32 +115,6 @@ end ##### A prescribed ocean... ##### -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} - architecture :: A - grid :: G - clock :: Clock{C} - velocities :: U - tracers :: T - timeseries :: F -end - -function PrescribedOcean(timeseries; - grid, - clock=Clock{Float64}(time = 0)) - - τx = Field{Face, Center, Nothing}(grid) - τy = Field{Center, Face, Nothing}(grid) - Jᵀ = Field{Center, Center, Nothing}(grid) - Jˢ = Field{Center, Center, Nothing}(grid) - - u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) - v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) - T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) - S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - - PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) -end - # ...with prescribed velocity and tracer fields version = ECCO4Monthly() dates = all_ECCO_dates(version)[1:24] @@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch) ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days) -##### -##### Need to extend a couple of methods -##### - -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) - tick!(model.clock, Δt) - time = Time(model.clock.time) - - possible_fts = merge(model.velocities, model.tracers) - time_series_tuple = extract_field_time_series(possible_fts) - - for fts in time_series_tuple - update_field_time_series!(fts, time) - end - - parent(model.velocities.u) .= parent(model.timeseries.u[time]) - parent(model.velocities.v) .= parent(model.timeseries.v[time]) - parent(model.tracers.T) .= parent(model.timeseries.T[time]) - parent(model.tracers.S) .= parent(model.timeseries.S[time]) - - return nothing -end - -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing - -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 - ##### ##### A prescribed atmosphere... ##### diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index 3ddaf5abd..e8251b4af 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -151,7 +151,9 @@ function regrid_bathymetry(target_grid; Nxn = length(λ_data) Nyn = length(φ_data) Nzn = 1 - + + @show φ₁_data, φ₂_data + @show λ₁_data, λ₂_data native_grid = LatitudeLongitudeGrid(arch, Float32; size = (Nxn, Nyn, Nzn), diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index dca368786..0dcd4ebac 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation +export ocean_simulation, PrescribedOcean using Oceananigans using Oceananigans.Units diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean.jl new file mode 100644 index 000000000..0a3e77f6f --- /dev/null +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -0,0 +1,71 @@ + +import Oceananigans.TimeSteppers: time_step!, update_state!, tick! +import Oceananigans.Models: timestepper, update_model_field_time_series! +import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity + +##### +##### A prescribed ocean... +##### + +struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} + architecture :: A + grid :: G + clock :: Clock{C} + velocities :: U + tracers :: T + timeseries :: F +end + +""" + PrescribedOcean(timeseries; grid, clock=Clock{Float64}(time = 0)) + +Create a `PrescribedOcean` model on a specified `grid` with a `clock` that evolves +according to the data passed in `timeseries`. + +`timeseries` _must_ be a `NamedTuple` containing `FieldTimeseries` for `u`, `v`, `T`, and `S`. +""" +function PrescribedOcean(timeseries; + grid, + clock=Clock{Float64}(time = 0)) + + τx = Field{Face, Center, Nothing}(grid) + τy = Field{Center, Face, Nothing}(grid) + Jᵀ = Field{Center, Center, Nothing}(grid) + Jˢ = Field{Center, Center, Nothing}(grid) + + u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) + v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) + T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) + S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) + + return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) +end + +##### +##### Need to extend a couple of methods +##### + +function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) + tick!(model.clock, Δt) + time = Time(model.clock.time) + + possible_fts = merge(model.velocities, model.tracers) + time_series_tuple = extract_field_time_series(possible_fts) + + for fts in time_series_tuple + update_field_time_series!(fts, time) + end + + parent(model.velocities.u) .= parent(model.timeseries.u[time]) + parent(model.velocities.v) .= parent(model.timeseries.v[time]) + parent(model.tracers.T) .= parent(model.timeseries.T[time]) + parent(model.tracers.S) .= parent(model.timeseries.S[time]) + + return nothing +end + +update_state!(::PrescribedOcean) = nothing +timestepper(::PrescribedOcean) = nothing + +reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 From b00046699da9a6c26ecb22449bfc5aba121a72fe Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 12:32:20 +0100 Subject: [PATCH 06/56] arctic simulation --- experiments/arctic_simulation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index d291e38b5..aa80e666d 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -36,10 +36,10 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) version = ECCO2Daily() dates = all_ECCO_dates(version)[1:200] -u = ECCOFieldTimeSeries(:u_velocity, version; dates, time_indices_in_memory=2) -v = ECCOFieldTimeSeries(:v_velocity, version; dates, time_indices_in_memory=2) -T = ECCOFieldTimeSeries(:temperature, version; dates, time_indices_in_memory=2) -S = ECCOFieldTimeSeries(:salinity, version; dates, time_indices_in_memory=2) +u = ECCOFieldTimeSeries(:u_velocity, version; grid, dates, time_indices_in_memory=10) +v = ECCOFieldTimeSeries(:v_velocity, version; grid, dates, time_indices_in_memory=10) +T = ECCOFieldTimeSeries(:temperature, version; grid, dates, time_indices_in_memory=10) +S = ECCOFieldTimeSeries(:salinity, version; grid, dates, time_indices_in_memory=10) ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model, Δt=10minutes, verbose=false) From 410722ccf33496a26369864e9f6f67ba2a65e4a1 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 12:40:09 +0100 Subject: [PATCH 07/56] works --- test/test_distributed_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_distributed_utils.jl b/test/test_distributed_utils.jl index 8c7906901..3910f06db 100644 --- a/test/test_distributed_utils.jl +++ b/test/test_distributed_utils.jl @@ -80,7 +80,7 @@ end @testset "Distributed Bathymetry interpolation" begin # We start by building a fake bathyemtry on rank 0 and save it to file @root begin - λ = -180:0.1:180 + λ = -179.5:0.1:179.5 φ = 0:0.1:50 Nλ = length(λ) From 5f343d68254c26f66b079f6554c95b9a407ff49a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 12:42:04 +0100 Subject: [PATCH 08/56] let's go --- src/ClimaOcean.jl | 2 +- src/OceanSimulations/OceanSimulations.jl | 1 + .../{prescribed_ocean.jl => prescribed_ocean_model.jl} | 5 +++++ 3 files changed, 7 insertions(+), 1 deletion(-) rename src/OceanSimulations/{prescribed_ocean.jl => prescribed_ocean_model.jl} (96%) diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 2205ede93..367dbd35c 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -73,9 +73,9 @@ end end include("DistributedUtils.jl") +include("OceanSeaIceModels/OceanSeaIceModels.jl") include("OceanSimulations/OceanSimulations.jl") include("SeaIceSimulations.jl") -include("OceanSeaIceModels/OceanSeaIceModels.jl") include("VerticalGrids.jl") include("InitialConditions/InitialConditions.jl") include("DataWrangling/DataWrangling.jl") diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 0dcd4ebac..02f517c1d 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -41,5 +41,6 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("ocean_simulation.jl") +include("prescribed_ocean_model.jl") end # module diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean_model.jl similarity index 96% rename from src/OceanSimulations/prescribed_ocean.jl rename to src/OceanSimulations/prescribed_ocean_model.jl index 0a3e77f6f..3433c4778 100644 --- a/src/OceanSimulations/prescribed_ocean.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -1,6 +1,11 @@ +using Oceananigans +using Oceananigans.BoundaryConditions +using Oceananigans.Grids: architecture + import Oceananigans.TimeSteppers: time_step!, update_state!, tick! import Oceananigans.Models: timestepper, update_model_field_time_series! + import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity ##### From 6dfa96d435f310a692b927ef54a29cf713656f06 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 12:52:35 +0100 Subject: [PATCH 09/56] continue --- experiments/arctic_simulation.jl | 4 ++-- src/Bathymetry.jl | 3 --- src/ClimaOcean.jl | 6 +++++- src/OceanSeaIceModels/OceanSeaIceModels.jl | 4 ++-- src/OceanSimulations/prescribed_ocean_model.jl | 4 +++- 5 files changed, 12 insertions(+), 9 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index aa80e666d..9e817e1e4 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -19,7 +19,7 @@ Nz = length(r_faces) - 1 grid = RotatedLatitudeLongitudeGrid(size = (Nx, Ny, Nz), latitude = (-45, 45), longitude = (-45, 45), - z = z_faces, + z = r_faces, north_pole = (180, 0), halo = (5, 5, 4), topology = (Bounded, Bounded, Bounded)) @@ -34,7 +34,7 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) # ...with prescribed velocity and tracer fields version = ECCO2Daily() -dates = all_ECCO_dates(version)[1:200] +dates = all_ECCO_dates(version)[1:30] u = ECCOFieldTimeSeries(:u_velocity, version; grid, dates, time_indices_in_memory=10) v = ECCOFieldTimeSeries(:v_velocity, version; grid, dates, time_indices_in_memory=10) diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index e8251b4af..9d4f1631f 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -152,9 +152,6 @@ function regrid_bathymetry(target_grid; Nyn = length(φ_data) Nzn = 1 - @show φ₁_data, φ₂_data - @show λ₁_data, λ₂_data - native_grid = LatitudeLongitudeGrid(arch, Float32; size = (Nxn, Nyn, Nzn), latitude = (φ₁_data, φ₂_data), diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 367dbd35c..9e3756576 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -43,6 +43,10 @@ using DataDeps using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries using Oceananigans.Grids: node +using SeawaterPolynomials.TEOS10 + +function reference_density end +function heat_capacity end const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries, GPUAdaptedFieldTimeSeries} @@ -73,9 +77,9 @@ end end include("DistributedUtils.jl") -include("OceanSeaIceModels/OceanSeaIceModels.jl") include("OceanSimulations/OceanSimulations.jl") include("SeaIceSimulations.jl") +include("OceanSeaIceModels/OceanSeaIceModels.jl") include("VerticalGrids.jl") include("InitialConditions/InitialConditions.jl") include("DataWrangling/DataWrangling.jl") diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 65d796588..e5fc73094 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -23,13 +23,13 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using ClimaOcean: stateindex +import ClimaOcean: reference_density, heat_capacity + using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll function downwelling_radiation end function freshwater_flux end -function reference_density end -function heat_capacity end const default_gravitational_acceleration = 9.80665 const default_freshwater_density = 1000 diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 3433c4778..efcf9624d 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -3,10 +3,12 @@ using Oceananigans using Oceananigans.BoundaryConditions using Oceananigans.Grids: architecture +using Oceananigans.Models: AbstractModel + import Oceananigans.TimeSteppers: time_step!, update_state!, tick! import Oceananigans.Models: timestepper, update_model_field_time_series! -import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity +import ClimaOcean: reference_density, heat_capacity ##### ##### A prescribed ocean... From 76bd832d8cc55fbb4693476fa3627f084f8b737b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 13:26:38 +0100 Subject: [PATCH 10/56] a default inpainting for ECCO --- experiments/arctic_simulation.jl | 13 +++++++++---- src/DataWrangling/ECCO/ECCO.jl | 3 ++- src/DataWrangling/ECCO/ECCO_restoring.jl | 2 +- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 9e817e1e4..a22ff89a7 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -36,10 +36,15 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) version = ECCO2Daily() dates = all_ECCO_dates(version)[1:30] -u = ECCOFieldTimeSeries(:u_velocity, version; grid, dates, time_indices_in_memory=10) -v = ECCOFieldTimeSeries(:v_velocity, version; grid, dates, time_indices_in_memory=10) -T = ECCOFieldTimeSeries(:temperature, version; grid, dates, time_indices_in_memory=10) -S = ECCOFieldTimeSeries(:salinity, version; grid, dates, time_indices_in_memory=10) +u_metadata = ECCOMetadata(:u_velocity; version, dates) +v_metadata = ECCOMetadata(:v_velocity; version, dates) +T_metadata = ECCOMetadata(:temperature; version, dates) +S_metadata = ECCOMetadata(:salinity; version, dates) + +u = ECCOFieldTimeSeries(u_metadata, grid; time_indices_in_memory=2) +v = ECCOFieldTimeSeries(v_metadata, grid; time_indices_in_memory=2) +T = ECCOFieldTimeSeries(T_metadata, grid; time_indices_in_memory=2) +S = ECCOFieldTimeSeries(S_metadata, grid; time_indices_in_memory=2) ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model, Δt=10minutes, verbose=false) diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl index be78f7d2f..120a0f753 100644 --- a/src/DataWrangling/ECCO/ECCO.jl +++ b/src/DataWrangling/ECCO/ECCO.jl @@ -32,7 +32,6 @@ end include("ECCO_metadata.jl") include("ECCO_mask.jl") -include("ECCO_restoring.jl") # Vertical coordinate const ECCO_z = [ @@ -273,5 +272,7 @@ function set!(field::Field, ECCO_metadata::ECCOMetadata; kw...) return field end +include("ECCO_restoring.jl") + end # Module diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl index a8b89ee74..78fcac463 100644 --- a/src/DataWrangling/ECCO/ECCO_restoring.jl +++ b/src/DataWrangling/ECCO/ECCO_restoring.jl @@ -147,7 +147,7 @@ end function ECCOFieldTimeSeries(metadata::ECCOMetadata, grid::AbstractGrid; time_indices_in_memory = 2, time_indexing = Cyclical(), - inpainting = nothing, + inpainting = default_inpainting(metadata), cache_inpainted_data = true) # Make sure all the required individual files are downloaded From b56a44c1efdc5b1c9a742e02d0e9871b18a9ba13 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 13:48:20 +0100 Subject: [PATCH 11/56] try it out --- experiments/arctic_simulation.jl | 72 ++++++++++++++++++++++++-------- 1 file changed, 55 insertions(+), 17 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index a22ff89a7..a02140c5a 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -29,25 +29,22 @@ bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) ##### -##### A Prescribed Ocean model +##### A Propgnostic Ocean model ##### -# ...with prescribed velocity and tracer fields -version = ECCO2Daily() -dates = all_ECCO_dates(version)[1:30] +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = Centered() -u_metadata = ECCOMetadata(:u_velocity; version, dates) -v_metadata = ECCOMetadata(:v_velocity; version, dates) -T_metadata = ECCOMetadata(:temperature; version, dates) -S_metadata = ECCOMetadata(:salinity; version, dates) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) -u = ECCOFieldTimeSeries(u_metadata, grid; time_indices_in_memory=2) -v = ECCOFieldTimeSeries(v_metadata, grid; time_indices_in_memory=2) -T = ECCOFieldTimeSeries(T_metadata, grid; time_indices_in_memory=2) -S = ECCOFieldTimeSeries(S_metadata, grid; time_indices_in_memory=2) +ocean = ocean_simulation(grid; + momentum_advection, + tracer_advection, + free_surface) -ocean_model = PrescribedOcean((; u, v, T, S); grid) -ocean = Simulation(ocean_model, Δt=10minutes, verbose=false) + +set!(ocean.model, T=ECCOMetadata(:temperature), + S=ECCOMetadata(:salinity)) ##### ##### A Prognostic Sea-ice model @@ -59,7 +56,7 @@ set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), ℵ=ECCOMetadata(:sea_ice_concentration)) ##### -##### Atmosphere model +##### A Prescribed Atmosphere model ##### atmosphere = JRA55PrescribedAtmosphere(; backend=JRA55NetCDFBackend(40)) @@ -80,7 +77,13 @@ Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu), filename = "sea_ice_quantities.jld2", - schedule = IterationInterval(12)) + schedule = IterationInterval(12), + overwrite_existing=true) + +sea_ice.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ), + filename = "averaged_sea_ice_quantities.jld2", + schedule = AveragedTimeInterval(1days), + overwrite_existing=true) wall_time = Ref(time_ns()) @@ -113,4 +116,39 @@ end # And add it as a callback to the simulation. add_callback!(arctic, progress, IterationInterval(10)) -run!(arctic) \ No newline at end of file +run!(arctic) + +##### +##### Comparison to ECCO Climatology +##### + +version = ECCO4Monthly() +dates = all_ECCO_dates(version)[1:12] + +h_metadata = ECCOMetadata(:sea_ice_thickness; version, dates) +ℵ_metadata = ECCOMetadata(:sea_ice_concentration; version, dates) + +# Montly averaged ECCO data +hE = ECCOFieldTimeSeries(h_metadata, grid; time_indices_in_memory=12) +ℵE = ECCOFieldTimeSeries(ℵ_metadata, grid; time_indices_in_memory=12) + +# Daily averaged Model output +hm = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "h") +ℵm = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "ℵ") + +# Montly average the model output +h̄m = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) +ℵ̄m = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) + +for (i, time) in enumerate(hE.times) + counter = 0 + for (j, t) in enumerate(hm.times) + if t ≤ time + h̄m[i] .+= hm[j] + ℵ̄m[i] .+= ℵm[j] + counter += 1 + end + end + h̄m[i] ./= counter + ℵ̄m[i] ./= counter +end \ No newline at end of file From d3e0c54a27959fb5f5998a8386d29540715d7d02 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 08:19:59 -0500 Subject: [PATCH 12/56] run it on the GPU --- experiments/arctic_simulation.jl | 19 +++++++++++-------- .../component_interfaces.jl | 4 ++-- .../time_step_ocean_sea_ice_model.jl | 2 +- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index a02140c5a..84ea21c0a 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -9,6 +9,9 @@ using ClimaOcean.ECCO using ClimaOcean.ECCO: all_ECCO_dates using Printf +using CUDA +CUDA.device!(1) + r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=3000) z_faces = MutableVerticalDiscretization(r_faces) @@ -16,13 +19,13 @@ Nx = 180 # longitudinal direction -> 250 points is about 1.5ᵒ resolution Ny = 180 # meridional direction -> same thing, 48 points is about 1.5ᵒ resolution Nz = length(r_faces) - 1 -grid = RotatedLatitudeLongitudeGrid(size = (Nx, Ny, Nz), - latitude = (-45, 45), - longitude = (-45, 45), - z = r_faces, - north_pole = (180, 0), - halo = (5, 5, 4), - topology = (Bounded, Bounded, Bounded)) +grid = RotatedLatitudeLongitudeGrid(GPU(), size = (Nx, Ny, Nz), + latitude = (-45, 45), + longitude = (-45, 45), + z = r_faces, + north_pole = (180, 0), + halo = (5, 5, 4), + topology = (Bounded, Bounded, Bounded)) bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) @@ -59,7 +62,7 @@ set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), ##### A Prescribed Atmosphere model ##### -atmosphere = JRA55PrescribedAtmosphere(; backend=JRA55NetCDFBackend(40)) +atmosphere = JRA55PrescribedAtmosphere(GPU(); backend=JRA55NetCDFBackend(40)) radiation = Radiation() ##### diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index f83e0c846..359976700 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -83,11 +83,11 @@ 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) - + kernel_parameters = interface_kernel_parameters(exchange_grid) launch!(arch, exchange_grid, kernel_parameters, _compute_fractional_indices!, frac_indices, exchange_grid, atmos_grid) diff --git a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl index b8dd29a3e..c1a9e512a 100644 --- a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl @@ -32,7 +32,7 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_ end sea_ice.Δt = Δt - time_step!(sea_ice) + thermodynamic_sea_ice_time_step!(coupled_model) end # TODO after ice time-step: From 9473d7d036faf50322fb51260c9de14b98f3cb77 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 08:31:14 -0500 Subject: [PATCH 13/56] churning on --- .../InterfaceComputations/component_interfaces.jl | 8 +++----- src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl | 2 +- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 359976700..d7591430d 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -81,11 +81,9 @@ 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] - frac_indices = OffsetArray(frac_indices_data, -1, -1, 0) + 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) frac_indices = on_architecture(arch, frac_indices) kernel_parameters = interface_kernel_parameters(exchange_grid) diff --git a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl index c1a9e512a..b8dd29a3e 100644 --- a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl @@ -32,7 +32,7 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_ end sea_ice.Δt = Δt - thermodynamic_sea_ice_time_step!(coupled_model) + time_step!(sea_ice) end # TODO after ice time-step: From 08b065e113b4751b4e1a47a6ffc8b35572f2cc5c Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 08:46:20 -0500 Subject: [PATCH 14/56] try it out like this --- experiments/arctic_simulation.jl | 2 +- src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 84ea21c0a..d686ac8ff 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -70,7 +70,7 @@ radiation = Radiation() ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=600, stop_time=365days) +arctic = Simulation(arctic, Δt=300, stop_time=365days) h = sea_ice.model.ice_thickness ℵ = sea_ice.model.ice_concentration diff --git a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl index b8dd29a3e..c1a9e512a 100644 --- a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl @@ -32,7 +32,7 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_ end sea_ice.Δt = Δt - time_step!(sea_ice) + thermodynamic_sea_ice_time_step!(coupled_model) end # TODO after ice time-step: From 0ebcf53a3113072d82ce3f21a2faf25666a747b8 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 08:53:48 -0500 Subject: [PATCH 15/56] remove this constant coefficient fluxes --- .../InterfaceComputations/component_interfaces.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index d7591430d..ca6953f3a 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -223,9 +223,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(), - atmosphere_sea_ice_flux_formulation = CoefficientBasedFluxes(drag_coefficient=2e-3, - heat_transfer_coefficient=1e-4, - vapor_flux_coefficient=1e-4), + atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(), atmosphere_ocean_interface_temperature = BulkTemperature(), atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean), atmosphere_sea_ice_interface_temperature = default_ai_temperature(sea_ice), From 2003a46bc64f4f833499781b092e8bc230458896 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 4 Mar 2025 15:31:06 +0100 Subject: [PATCH 16/56] add fluxes --- experiments/arctic_simulation.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index a02140c5a..acd9a7b84 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -35,7 +35,7 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) momentum_advection = WENOVectorInvariant(order=5) tracer_advection = Centered() -free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.8) ocean = ocean_simulation(grid; momentum_advection, @@ -67,20 +67,26 @@ radiation = Radiation() ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=600, stop_time=365days) +arctic = Simulation(arctic, Δt=8minutes, stop_time=365days) h = sea_ice.model.ice_thickness ℵ = sea_ice.model.ice_concentration Gh = sea_ice.model.timestepper.Gⁿ.h Gℵ = sea_ice.model.timestepper.Gⁿ.ℵ Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature - -sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu), +Qˡ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat +Qˢ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat +Qⁱ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat +Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat +Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat +Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat + +sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "sea_ice_quantities.jld2", schedule = IterationInterval(12), overwrite_existing=true) -sea_ice.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ), +sea_ice.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "averaged_sea_ice_quantities.jld2", schedule = AveragedTimeInterval(1days), overwrite_existing=true) From 983565a67c9694c4441d00255222bb4175935e6a Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 09:33:45 -0500 Subject: [PATCH 17/56] a very diffusive ocean --- experiments/arctic_simulation.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 80f620116..5ef97e582 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -35,8 +35,9 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) ##### A Propgnostic Ocean model ##### -momentum_advection = WENOVectorInvariant(order=5) -tracer_advection = Centered() +# A very diffusive ocean +momentum_advection = WENOVectorInvariant(order=3) +tracer_advection = WENO(order=3) free_surface = SplitExplicitFreeSurface(grid; cfl=0.8) @@ -70,7 +71,7 @@ radiation = Radiation() ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=8minutes, stop_time=365days) +arctic = Simulation(arctic, Δt=20minutes, stop_time=365days) h = sea_ice.model.ice_thickness ℵ = sea_ice.model.ice_concentration From 30502c78f69348d9957389cf634a87f3f8bfbdc6 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 09:36:23 -0500 Subject: [PATCH 18/56] just push it with the diffusivity --- experiments/arctic_simulation.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 5ef97e582..bf8349a47 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -36,15 +36,17 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) ##### # A very diffusive ocean -momentum_advection = WENOVectorInvariant(order=3) -tracer_advection = WENO(order=3) +momentum_advection = VectorInvariant() +tracer_advection = Centered(order=3) free_surface = SplitExplicitFreeSurface(grid; cfl=0.8) +closure = (ClimaOcean.OceanSimulations.default_ocean_closure(), HorizontalScalarDiffusivity(κ=1000, ν=10000)) ocean = ocean_simulation(grid; momentum_advection, tracer_advection, - free_surface) + free_surface, + closure) set!(ocean.model, T=ECCOMetadata(:temperature), From a96e1055a87ec1257ba16a63b83c8720e1db52ad Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 09:53:51 -0500 Subject: [PATCH 19/56] unfortunately it crashes about 5 minutes --- experiments/arctic_simulation.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index bf8349a47..40ae0a104 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -12,7 +12,7 @@ using Printf using CUDA CUDA.device!(1) -r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=3000) +r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=2000) z_faces = MutableVerticalDiscretization(r_faces) Nx = 180 # longitudinal direction -> 250 points is about 1.5ᵒ resolution @@ -37,10 +37,10 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) # A very diffusive ocean momentum_advection = VectorInvariant() -tracer_advection = Centered(order=3) +tracer_advection = WENO(order=3) -free_surface = SplitExplicitFreeSurface(grid; cfl=0.8) -closure = (ClimaOcean.OceanSimulations.default_ocean_closure(), HorizontalScalarDiffusivity(κ=1000, ν=10000)) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) +closure = (ClimaOcean.OceanSimulations.default_ocean_closure(), HorizontalScalarDiffusivity(κ=1000, ν=3000)) ocean = ocean_simulation(grid; momentum_advection, @@ -73,12 +73,15 @@ radiation = Radiation() ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=20minutes, stop_time=365days) +arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) -h = sea_ice.model.ice_thickness -ℵ = sea_ice.model.ice_concentration +# Sea-ice variables +h = sea_ice.model.ice_thickness +ℵ = sea_ice.model.ice_concentration Gh = sea_ice.model.timestepper.Gⁿ.h Gℵ = sea_ice.model.timestepper.Gⁿ.ℵ + +# Fluxes Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature Qˡ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat Qˢ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat @@ -87,6 +90,7 @@ Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat +# Ou sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "sea_ice_quantities.jld2", schedule = IterationInterval(12), From 5b897af681b28dac34a78f1160d6659b074db618 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 10:14:56 -0500 Subject: [PATCH 20/56] try with 8mins --- experiments/arctic_simulation.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 40ae0a104..d002d8975 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -36,11 +36,11 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) ##### # A very diffusive ocean -momentum_advection = VectorInvariant() +momentum_advection = WENOVectorInvariant(order=3) tracer_advection = WENO(order=3) free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) -closure = (ClimaOcean.OceanSimulations.default_ocean_closure(), HorizontalScalarDiffusivity(κ=1000, ν=3000)) +closure = ClimaOcean.OceanSimulations.default_ocean_closure() ocean = ocean_simulation(grid; momentum_advection, @@ -48,7 +48,6 @@ ocean = ocean_simulation(grid; free_surface, closure) - set!(ocean.model, T=ECCOMetadata(:temperature), S=ECCOMetadata(:salinity)) @@ -73,7 +72,7 @@ radiation = Radiation() ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) +arctic = Simulation(arctic, Δt=8minutes, stop_time=365days) # Sea-ice variables h = sea_ice.model.ice_thickness @@ -90,7 +89,7 @@ Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat -# Ou +# Output writers sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "sea_ice_quantities.jld2", schedule = IterationInterval(12), From a512a0cf685a11aa270acea2ec6ffc9ceddc26e7 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Tue, 4 Mar 2025 13:25:59 -0500 Subject: [PATCH 21/56] to change evenutally --- examples/idealized_single_column_simulation.jl | 6 +++--- experiments/arctic_simulation.jl | 13 +++++++++---- .../component_interfaces.jl | 16 +++++++++------- .../sea_ice_ocean_fluxes.jl | 2 +- src/OceanSeaIceModels/ocean_sea_ice_model.jl | 6 +++--- .../time_step_ocean_sea_ice_model.jl | 4 ++-- src/SeaIceSimulations.jl | 12 ++++++------ test/test_simulations.jl | 2 +- 8 files changed, 34 insertions(+), 27 deletions(-) diff --git a/examples/idealized_single_column_simulation.jl b/examples/idealized_single_column_simulation.jl index d1cab50ce..41efed5ae 100644 --- a/examples/idealized_single_column_simulation.jl +++ b/examples/idealized_single_column_simulation.jl @@ -42,7 +42,7 @@ set!(ocean.model, T=Tᵢ, S=S₀) sea_ice_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat)) top_sea_ice_temperature = Field{Center, Center, Nothing}(sea_ice_grid) top_heat_boundary_condition = PrescribedTemperature(top_sea_ice_temperature) -ice_thermodynamics = SlabSeaIceThermodynamics(sea_ice_grid; top_heat_boundary_condition) +thermodynamics = SlabSeaIceThermodynamics(sea_ice_grid; top_heat_boundary_condition) top_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid) bottom_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid) @@ -50,7 +50,7 @@ bottom_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid) sea_ice_model = SeaIceModel(sea_ice_grid; top_heat_flux = top_sea_ice_heat_flux, bottom_heat_flux = bottom_sea_ice_heat_flux, - ice_thermodynamics) + thermodynamics) set!(sea_ice_model.ice_concentration, 0.9) set!(sea_ice_model.ice_thickness, 0.1) @@ -76,7 +76,7 @@ Qi = sea_ice_model.external_heat_fluxes.top fluxes = (; Qo, Qi, Qib, Qv, Qc, τx, τy, Fv) ocean_outputs = merge(ocean.model.velocities, ocean.model.tracers) sea_ice_outputs = (h = sea_ice.model.ice_thickness, - Ti = sea_ice.model.ice_thermodynamics.top_surface_temperature, + Ti = sea_ice.model.thermodynamics.top_surface_temperature, ℵ = sea_ice.model.ice_concentration) outputs = merge(ocean_outputs, sea_ice_outputs, fluxes) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index d002d8975..969e1e635 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -7,6 +7,7 @@ using Oceananigans.OrthogonalSphericalShellGrids using ClimaOcean.OceanSimulations using ClimaOcean.ECCO using ClimaOcean.ECCO: all_ECCO_dates +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium using Printf using CUDA @@ -55,12 +56,16 @@ set!(ocean.model, T=ECCOMetadata(:temperature), ##### A Prognostic Sea-ice model ##### -sea_ice = sea_ice_simulation(grid; dynamics=nothing, advection=nothing) +# Remember to pass the SSS as a bottom bc to the sea ice! +SSS = view(ocean.model.tracers.S, :, :, grid.Nz) +bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS) + +sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics=nothing, advection=nothing) set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), ℵ=ECCOMetadata(:sea_ice_concentration)) -##### +##### ##### A Prescribed Atmosphere model ##### @@ -90,12 +95,12 @@ Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat # Output writers -sea_ice.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), +arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "sea_ice_quantities.jld2", schedule = IterationInterval(12), overwrite_existing=true) -sea_ice.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), +arctic.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "averaged_sea_ice_quantities.jld2", schedule = AveragedTimeInterval(1days), overwrite_existing=true) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index ca6953f3a..14e251e56 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -81,8 +81,8 @@ function StateExchanger(ocean::Simulation, atmosphere) Nx, Ny, Nz = size(exchange_grid) # 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] + FT = eltype(atmos_grid) + frac_indices = [FractionalIndices(one(FT), one(FT), nothing) for i=1:Nx+2, j=1:Ny+2, k=1:1] frac_indices = OffsetArray(frac_indices, -1, -1, 0) frac_indices = on_architecture(arch, frac_indices) @@ -169,7 +169,7 @@ function atmosphere_sea_ice_interface(sea_ice::SeaIceSimulation, specific_humidity_formulation, temperature_formulation) - interface_temperature = sea_ice.model.ice_thermodynamics.top_surface_temperature + interface_temperature = sea_ice.model.thermodynamics.top_surface_temperature return AtmosphereInterface(fluxes, ai_flux_formulation, interface_temperature, properties) end @@ -204,7 +204,7 @@ end default_ai_temperature(sea_ice) = nothing function default_ai_temperature(sea_ice::SeaIceSimulation) - conductive_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux.parameters.flux + conductive_flux = sea_ice.model.thermodynamics.internal_heat_flux.parameters.flux return SkinTemperature(conductive_flux) end @@ -223,7 +223,9 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(), - atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(), + atmosphere_sea_ice_flux_formulation = CoefficientBasedFluxes(drag_coefficient=2e-3, + heat_transfer_coefficient=1e-4, + vapor_flux_coefficient=1e-4), atmosphere_ocean_interface_temperature = BulkTemperature(), atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean), atmosphere_sea_ice_interface_temperature = default_ai_temperature(sea_ice), @@ -267,7 +269,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; sea_ice_properties = (reference_density = sea_ice_reference_density, heat_capacity = sea_ice_heat_capacity, freshwater_density = freshwater_density, - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus, + liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus, temperature_units = sea_ice_temperature_units) net_top_sea_ice_fluxes = (; heat=sea_ice.model.external_heat_fluxes.top) @@ -313,7 +315,7 @@ function sea_ice_similarity_theory(sea_ice::SeaIceSimulation) # SkinTemperature. Also we need to pass the sea ice internal flux # The thickness and salinity need to be passed as well, # but the can be passed as state variables once we refactor the `StateValues` struct. - internal_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux + internal_flux = sea_ice.model.thermodynamics.internal_heat_flux interface_temperature_type = SkinTemperature(internal_flux) return SimilarityTheoryFluxes(; interface_temperature_type) end diff --git a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index 4d6d32577..1fbfc27d9 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -21,7 +21,7 @@ function compute_sea_ice_ocean_latent_heat_flux!(coupled_model) ℵᵢ = sea_ice.model.ice_concentration ocean_properties = coupled_model.interfaces.ocean_properties - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus grid = ocean.model.grid arch = architecture(grid) diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 7d829f53c..f66c30e6b 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -77,11 +77,11 @@ heat_capacity(unsupported) = reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation) reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) reference_density(eos::TEOS10EquationOfState) = eos.reference_density -reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density +reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.thermodynamics.phase_transitions.ice_density heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation) heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) -heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity +heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.thermodynamics.phase_transitions.ice_heat_capacity # Does not really matter if there is no model reference_density(::Nothing) = 0 @@ -180,7 +180,7 @@ function above_freezing_ocean_temperature!(ocean, sea_ice::SeaIceSimulation) T = ocean.model.tracers.T S = ocean.model.tracers.S ℵ = sea_ice.model.ice_concentration - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus grid = ocean.model.grid arch = architecture(grid) diff --git a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl index c1a9e512a..5c091052d 100644 --- a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl @@ -32,7 +32,7 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_ end sea_ice.Δt = Δt - thermodynamic_sea_ice_time_step!(coupled_model) + time_step!(sea_ice) end # TODO after ice time-step: @@ -77,7 +77,7 @@ function thermodynamic_sea_ice_time_step!(coupled_model) grid = coupled_model.sea_ice.model.grid arch = architecture(grid) clock = model.clock - thermodynamics = model.ice_thermodynamics + thermodynamics = model.thermodynamics ice_thickness = model.ice_thickness ice_concentration = model.ice_concentration ice_consolidation_thickness = model.ice_consolidation_thickness diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 834be7db1..f4f8e18ef 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -41,11 +41,11 @@ function sea_ice_simulation(grid; top_surface_temperature = Field{Center, Center, Nothing}(grid) top_heat_boundary_condition = PrescribedTemperature(top_surface_temperature) - ice_thermodynamics = SlabSeaIceThermodynamics(grid; - internal_heat_flux, - phase_transitions, - top_heat_boundary_condition, - bottom_heat_boundary_condition) + thermodynamics = SlabSeaIceThermodynamics(grid; + internal_heat_flux, + phase_transitions, + top_heat_boundary_condition, + bottom_heat_boundary_condition) bottom_heat_flux = Field{Center, Center, Nothing}(grid) top_heat_flux = Field{Center, Center, Nothing}(grid) @@ -60,7 +60,7 @@ function sea_ice_simulation(grid; tracers, ice_consolidation_thickness, # top_momentum_stress, - ice_thermodynamics, + thermodynamics, dynamics, bottom_heat_flux, top_heat_flux) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 96868eb5c..0965a58f9 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -48,7 +48,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature sea_ice_grid = ImmersedBoundaryGrid(sea_ice_grid, GridFittedBottom(bottom_height)) sea_ice = sea_ice_simulation(sea_ice_grid) - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus # Set the ocean temperature and salinity set!(ocean.model, T=temperature_metadata[1], S=salinity_metadata[1]) From 72abc67d952909b8acf7cce5993d0bb8b3952860 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 5 Mar 2025 11:14:09 -0500 Subject: [PATCH 22/56] 5 minutes --- experiments/arctic_simulation.jl | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 969e1e635..135ccdc0e 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -77,7 +77,7 @@ radiation = Radiation() ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=8minutes, stop_time=365days) +arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) # Sea-ice variables h = sea_ice.model.ice_thickness @@ -153,22 +153,23 @@ hE = ECCOFieldTimeSeries(h_metadata, grid; time_indices_in_memory=12) ℵE = ECCOFieldTimeSeries(ℵ_metadata, grid; time_indices_in_memory=12) # Daily averaged Model output -hm = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "h") -ℵm = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "ℵ") +h = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "h") +ℵ = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "ℵ") # Montly average the model output -h̄m = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) -ℵ̄m = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) +hm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) +ℵm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) -for (i, time) in enumerate(hE.times) +for (i, time) in enumerate(hm.times) counter = 0 - for (j, t) in enumerate(hm.times) + for (j, t) in enumerate(h.times) if t ≤ time - h̄m[i] .+= hm[j] - ℵ̄m[i] .+= ℵm[j] + hm[i] .+= h[j] + ℵm[i] .+= ℵ[j] counter += 1 end end - h̄m[i] ./= counter - ℵ̄m[i] ./= counter -end \ No newline at end of file + @show counter + hm[i] ./= counter + ℵm[i] ./= counter +end From cd7cec42ed9791d97daac649bb817254ae65cfc7 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 5 Mar 2025 11:49:22 -0500 Subject: [PATCH 23/56] add a momentum equation --- experiments/arctic_simulation.jl | 20 ++++++++++++++++++- .../assemble_net_fluxes.jl | 13 ++++++++---- .../component_interfaces.jl | 12 ++++++++++- src/SeaIceSimulations.jl | 4 ---- 4 files changed, 39 insertions(+), 10 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 135ccdc0e..6745293f8 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -56,11 +56,29 @@ set!(ocean.model, T=ECCOMetadata(:temperature), ##### A Prognostic Sea-ice model ##### +using ClimaSeaIce.SeaIceMomentumEquations +using ClimaSeaIce.Rheologies + # Remember to pass the SSS as a bottom bc to the sea ice! SSS = view(ocean.model.tracers.S, :, :, grid.Nz) bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS) -sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics=nothing, advection=nothing) +SSU = view(ocean.model.velocities.u, :, :, grid.Nz) +SSV = view(ocean.model.velocities.u, :, :, grid.Nz) + +τo = SemiImplicitStress(uₑ=SSU, vₑ=SSV) +τua = Field{Face, Center, Nothing}(grid) +τva = Field{Center, Face, Nothing}(grid) + +dynamics = SeaIceMomentumEquation(grid; + coriolis = ocean.model.coriolis, + top_momentum_stress = (u=τua, v=τva), + bottom_momentum_stress = τo, + ocean_velocities = (u=SSU, v=SSV), + rhelogy = ElastoViscoPlasticRheology(), + solver = SplitExplicitSolver(120)) + +sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), ℵ=ECCOMetadata(:sea_ice_concentration)) diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 07ba1a901..37bbb29c9 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -154,7 +154,7 @@ function compute_net_sea_ice_fluxes!(coupled_model) arch = architecture(grid) clock = coupled_model.clock - top_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_top.heat + top_fluxes = coupled_model.interfaces.net_fluxes.sea_ice_top bottom_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_bottom.heat sea_ice_ocean_fluxes = coupled_model.interfaces.sea_ice_ocean_interface.fluxes atmosphere_sea_ice_fluxes = coupled_model.interfaces.atmosphere_sea_ice_interface.fluxes @@ -177,7 +177,7 @@ function compute_net_sea_ice_fluxes!(coupled_model) launch!(arch, grid, kernel_parameters, _assemble_net_sea_ice_fluxes!, - top_heat_flux, + top_fluxes, bottom_heat_flux, grid, clock, @@ -192,7 +192,7 @@ function compute_net_sea_ice_fluxes!(coupled_model) return nothing end -@kernel function _assemble_net_sea_ice_fluxes!(top_heat_flux, +@kernel function _assemble_net_sea_ice_fluxes!(top_fluxes, bottom_heat_flux, grid, clock, @@ -220,6 +220,9 @@ end Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1] # interfacial heat flux end + ρτx = atmosphere_sea_ice_fluxes.x_momentum # zonal momentum flux + ρτy = atmosphere_sea_ice_fluxes.y_momentum # meridional momentum flux + # Compute radiation fluxes σ = atmos_sea_ice_properties.radiation.σ α = stateindex(atmos_sea_ice_properties.radiation.α, i, j, kᴺ, grid, time) @@ -233,6 +236,8 @@ end # Mask fluxes over land for convenience inactive = inactive_node(i, j, kᴺ, grid, c, c, c) - @inbounds top_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQt) + @inbounds top_fluxes.heat[i, j, 1] = ifelse(inactive, zero(grid), ΣQt) + @inbounds top_fluxes.u[i, j, 1] = ifelse(inactive, zero(grid), ℑxᶠᵃᵃ(i, j, 1, grid, ρτx)) + @inbounds top_fluxes.v[i, j, 1] = ifelse(inactive, zero(grid), ℑyᵃᶠᵃ(i, j, 1, grid, ρτy)) @inbounds bottom_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQb) end diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 14e251e56..7b4abeed3 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -272,7 +272,17 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus, temperature_units = sea_ice_temperature_units) - net_top_sea_ice_fluxes = (; heat=sea_ice.model.external_heat_fluxes.top) + net_momentum_fluxes = if sea_ice.model.dynamics isa Nothing + u = Field{Face, Center, Nothing}(sea_ice.model.grid) + v = Field{Center, Face, Nothing}(sea_ice.model.grid) + (; u, v) + else + u = sea_ice.model.dynamics.external_momentum_fluxes.top.u + v = sea_ice.model.dynamics.external_momentum_fluxes.top.v + (u, v) + end + + net_top_sea_ice_fluxes = merge((; heat=sea_ice.model.external_heat_fluxes.top), net_momentum_fluxes) net_bottom_sea_ice_fluxes = (; heat=sea_ice.model.external_heat_fluxes.bottom) else sea_ice_properties = nothing diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index f4f8e18ef..76e78231d 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -50,16 +50,12 @@ function sea_ice_simulation(grid; bottom_heat_flux = Field{Center, Center, Nothing}(grid) top_heat_flux = Field{Center, Center, Nothing}(grid) - # top_momentum_stress = (u = Field{Face, Center, Nothing}(grid), - # v = Field{Center, Face, Nothing}(grid)) - # Build the sea ice model sea_ice_model = SeaIceModel(grid; ice_salinity, advection, tracers, ice_consolidation_thickness, - # top_momentum_stress, thermodynamics, dynamics, bottom_heat_flux, From 394dce65af3b6980db2aeebc9f82491cc9ac849b Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 5 Mar 2025 11:51:00 -0500 Subject: [PATCH 24/56] with momentum --- experiments/arctic_simulation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 6745293f8..793d66c28 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -75,7 +75,7 @@ dynamics = SeaIceMomentumEquation(grid; top_momentum_stress = (u=τua, v=τva), bottom_momentum_stress = τo, ocean_velocities = (u=SSU, v=SSV), - rhelogy = ElastoViscoPlasticRheology(), + rheology = ElastoViscoPlasticRheology(), solver = SplitExplicitSolver(120)) sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) From 4fd4a629bcea501deee120c880344741559b5b5e Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Thu, 6 Mar 2025 07:15:37 -0500 Subject: [PATCH 25/56] add stuff --- experiments/arctic_simulation.jl | 14 ++++++++------ .../InterfaceComputations/component_interfaces.jl | 6 +++--- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 793d66c28..c75638df7 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -98,10 +98,10 @@ arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) # Sea-ice variables -h = sea_ice.model.ice_thickness -ℵ = sea_ice.model.ice_concentration -Gh = sea_ice.model.timestepper.Gⁿ.h -Gℵ = sea_ice.model.timestepper.Gⁿ.ℵ +h = sea_ice.model.ice_thickness +ℵ = sea_ice.model.ice_concentration +u = sea_ice.model.velocities.u +v = sea_ice.model.velocities.v # Fluxes Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature @@ -111,14 +111,16 @@ Qⁱ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat +τx = arctic.model.interfaces.net_fluxes.sea_ice_top.u +τy = arctic.model.interfaces.net_fluxes.sea_ice_top.v # Output writers -arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), +arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, u, v, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, τx, τy), filename = "sea_ice_quantities.jld2", schedule = IterationInterval(12), overwrite_existing=true) -arctic.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), +arctic.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, u, v, τx, τy), filename = "averaged_sea_ice_quantities.jld2", schedule = AveragedTimeInterval(1days), overwrite_existing=true) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 7b4abeed3..c9713c657 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -277,9 +277,9 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; v = Field{Center, Face, Nothing}(sea_ice.model.grid) (; u, v) else - u = sea_ice.model.dynamics.external_momentum_fluxes.top.u - v = sea_ice.model.dynamics.external_momentum_fluxes.top.v - (u, v) + u = sea_ice.model.dynamics.external_momentum_stresses.top.u + v = sea_ice.model.dynamics.external_momentum_stresses.top.v + (; u, v) end net_top_sea_ice_fluxes = merge((; heat=sea_ice.model.external_heat_fluxes.top), net_momentum_fluxes) From c1a406b62c64c9f484a4c7b360c6af8da322bfa4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 16:49:04 +0100 Subject: [PATCH 26/56] vestigial code --- .../idealized_single_column_simulation.jl | 6 +- .../flux_climatology/flux_climatology.jl | 57 ++++++++++++++++++- src/Bathymetry.jl | 1 - src/ClimaOcean.jl | 4 -- .../sea_ice_ocean_fluxes.jl | 2 +- 5 files changed, 60 insertions(+), 10 deletions(-) diff --git a/examples/idealized_single_column_simulation.jl b/examples/idealized_single_column_simulation.jl index 41efed5ae..d1cab50ce 100644 --- a/examples/idealized_single_column_simulation.jl +++ b/examples/idealized_single_column_simulation.jl @@ -42,7 +42,7 @@ set!(ocean.model, T=Tᵢ, S=S₀) sea_ice_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat)) top_sea_ice_temperature = Field{Center, Center, Nothing}(sea_ice_grid) top_heat_boundary_condition = PrescribedTemperature(top_sea_ice_temperature) -thermodynamics = SlabSeaIceThermodynamics(sea_ice_grid; top_heat_boundary_condition) +ice_thermodynamics = SlabSeaIceThermodynamics(sea_ice_grid; top_heat_boundary_condition) top_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid) bottom_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid) @@ -50,7 +50,7 @@ bottom_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid) sea_ice_model = SeaIceModel(sea_ice_grid; top_heat_flux = top_sea_ice_heat_flux, bottom_heat_flux = bottom_sea_ice_heat_flux, - thermodynamics) + ice_thermodynamics) set!(sea_ice_model.ice_concentration, 0.9) set!(sea_ice_model.ice_thickness, 0.1) @@ -76,7 +76,7 @@ Qi = sea_ice_model.external_heat_fluxes.top fluxes = (; Qo, Qi, Qib, Qv, Qc, τx, τy, Fv) ocean_outputs = merge(ocean.model.velocities, ocean.model.tracers) sea_ice_outputs = (h = sea_ice.model.ice_thickness, - Ti = sea_ice.model.thermodynamics.top_surface_temperature, + Ti = sea_ice.model.ice_thermodynamics.top_surface_temperature, ℵ = sea_ice.model.ice_concentration) outputs = merge(ocean_outputs, sea_ice_outputs, fluxes) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index aff118f56..94a1a951e 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -13,7 +13,7 @@ using Adapt using Printf using KernelAbstractions: @index, @kernel -import Oceananigans.TimeSteppers: time_step!, update_state! +import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! import Oceananigans.Models: timestepper, update_model_field_time_series! import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity @@ -115,6 +115,32 @@ end ##### A prescribed ocean... ##### +struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} + architecture :: A + grid :: G + clock :: Clock{C} + velocities :: U + tracers :: T + timeseries :: F +end + +function PrescribedOcean(timeseries; + grid, + clock=Clock{Float64}(time = 0)) + + τx = Field{Face, Center, Nothing}(grid) + τy = Field{Center, Face, Nothing}(grid) + Jᵀ = Field{Center, Center, Nothing}(grid) + Jˢ = Field{Center, Center, Nothing}(grid) + + u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) + v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) + T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) + S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) + + PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) +end + # ...with prescribed velocity and tracer fields version = ECCO4Monthly() dates = all_ECCO_dates(version)[1:24] @@ -131,6 +157,35 @@ grid = ECCO.ECCO_immersed_grid(arch) ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days) +##### +##### Need to extend a couple of methods +##### + +function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) + tick!(model.clock, Δt) + time = Time(model.clock.time) + + possible_fts = merge(model.velocities, model.tracers) + time_series_tuple = extract_field_time_series(possible_fts) + + for fts in time_series_tuple + update_field_time_series!(fts, time) + end + + parent(model.velocities.u) .= parent(model.timeseries.u[time]) + parent(model.velocities.v) .= parent(model.timeseries.v[time]) + parent(model.tracers.T) .= parent(model.timeseries.T[time]) + parent(model.tracers.S) .= parent(model.timeseries.S[time]) + + return nothing +end + +update_state!(::PrescribedOcean) = nothing +timestepper(::PrescribedOcean) = nothing + +reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 + ##### ##### A prescribed atmosphere... ##### diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index 4a1ec8de0..753d6dc94 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -143,7 +143,6 @@ function regrid_bathymetry(target_grid; close(dataset) # Diagnose target grid information - arch = architecture(target_grid) arch = architecture(target_grid) λ_data = λ_data |> Array{BigFloat} φ_data = φ_data |> Array{BigFloat} diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index bb1c51d68..670149b21 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -44,10 +44,6 @@ using DataDeps using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries using Oceananigans.Grids: node -using SeawaterPolynomials.TEOS10 - -function reference_density end -function heat_capacity end const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries, GPUAdaptedFieldTimeSeries} diff --git a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index 914c72913..72769e130 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -21,7 +21,7 @@ function compute_sea_ice_ocean_latent_heat_flux!(coupled_model) ℵᵢ = sea_ice.model.ice_concentration ocean_properties = coupled_model.interfaces.ocean_properties - liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus grid = ocean.model.grid arch = architecture(grid) From d7bfc840fea91a1b6df9da54c8a4e4af3e72eab2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 16:50:05 +0100 Subject: [PATCH 27/56] vestigial code --- src/OceanSeaIceModels/OceanSeaIceModels.jl | 4 +- src/OceanSeaIceModels/ocean_sea_ice_model.jl | 6 +- src/OceanSimulations/OceanSimulations.jl | 3 +- .../prescribed_ocean_model.jl | 78 ------------------- src/SeaIceSimulations.jl | 16 ++-- test/test_simulations.jl | 2 +- 6 files changed, 17 insertions(+), 92 deletions(-) delete mode 100644 src/OceanSimulations/prescribed_ocean_model.jl diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 1e5e2d2dc..8041f497d 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -23,13 +23,13 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using ClimaOcean: stateindex -import ClimaOcean: reference_density, heat_capacity - using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll function downwelling_radiation end function freshwater_flux end +function reference_density end +function heat_capacity end const default_gravitational_acceleration = 9.80665 const default_freshwater_density = 1000 diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 5422336bd..d2344dabb 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -82,11 +82,11 @@ heat_capacity(unsupported) = reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation) reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) reference_density(eos::TEOS10EquationOfState) = eos.reference_density -reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.thermodynamics.phase_transitions.ice_density +reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation) heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) -heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.thermodynamics.phase_transitions.ice_heat_capacity +heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity # Does not really matter if there is no model reference_density(::Nothing) = 0 @@ -183,7 +183,7 @@ function above_freezing_ocean_temperature!(ocean, sea_ice::SeaIceSimulation) T = ocean.model.tracers.T S = ocean.model.tracers.S ℵ = sea_ice.model.ice_concentration - liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus grid = ocean.model.grid arch = architecture(grid) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 87aa5dd62..32423c690 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation, PrescribedOcean +export ocean_simulation using Oceananigans using Oceananigans.Units @@ -41,6 +41,5 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("ocean_simulation.jl") -include("prescribed_ocean_model.jl") end # module diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl deleted file mode 100644 index efcf9624d..000000000 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ /dev/null @@ -1,78 +0,0 @@ - -using Oceananigans -using Oceananigans.BoundaryConditions -using Oceananigans.Grids: architecture - -using Oceananigans.Models: AbstractModel - -import Oceananigans.TimeSteppers: time_step!, update_state!, tick! -import Oceananigans.Models: timestepper, update_model_field_time_series! - -import ClimaOcean: reference_density, heat_capacity - -##### -##### A prescribed ocean... -##### - -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} - architecture :: A - grid :: G - clock :: Clock{C} - velocities :: U - tracers :: T - timeseries :: F -end - -""" - PrescribedOcean(timeseries; grid, clock=Clock{Float64}(time = 0)) - -Create a `PrescribedOcean` model on a specified `grid` with a `clock` that evolves -according to the data passed in `timeseries`. - -`timeseries` _must_ be a `NamedTuple` containing `FieldTimeseries` for `u`, `v`, `T`, and `S`. -""" -function PrescribedOcean(timeseries; - grid, - clock=Clock{Float64}(time = 0)) - - τx = Field{Face, Center, Nothing}(grid) - τy = Field{Center, Face, Nothing}(grid) - Jᵀ = Field{Center, Center, Nothing}(grid) - Jˢ = Field{Center, Center, Nothing}(grid) - - u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) - v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) - T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) - S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - - return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) -end - -##### -##### Need to extend a couple of methods -##### - -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) - tick!(model.clock, Δt) - time = Time(model.clock.time) - - possible_fts = merge(model.velocities, model.tracers) - time_series_tuple = extract_field_time_series(possible_fts) - - for fts in time_series_tuple - update_field_time_series!(fts, time) - end - - parent(model.velocities.u) .= parent(model.timeseries.u[time]) - parent(model.velocities.v) .= parent(model.timeseries.v[time]) - parent(model.tracers.T) .= parent(model.timeseries.T[time]) - parent(model.tracers.S) .= parent(model.timeseries.S[time]) - - return nothing -end - -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing - -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 23410012a..2b7dc17b8 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -41,22 +41,26 @@ function sea_ice_simulation(grid; top_surface_temperature = Field{Center, Center, Nothing}(grid) top_heat_boundary_condition = PrescribedTemperature(top_surface_temperature) - thermodynamics = SlabSeaIceThermodynamics(grid; - internal_heat_flux, - phase_transitions, - top_heat_boundary_condition, - bottom_heat_boundary_condition) + ice_thermodynamics = SlabSeaIceThermodynamics(grid; + internal_heat_flux, + phase_transitions, + top_heat_boundary_condition, + bottom_heat_boundary_condition) bottom_heat_flux = Field{Center, Center, Nothing}(grid) top_heat_flux = Field{Center, Center, Nothing}(grid) + # top_momentum_stress = (u = Field{Face, Center, Nothing}(grid), + # v = Field{Center, Face, Nothing}(grid)) + # Build the sea ice model sea_ice_model = SeaIceModel(grid; ice_salinity, advection, tracers, ice_consolidation_thickness, - thermodynamics, + # top_momentum_stress, + ice_thermodynamics, dynamics, bottom_heat_flux, top_heat_flux) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 52e241051..a892b427b 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -50,7 +50,7 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature sea_ice_grid = ImmersedBoundaryGrid(sea_ice_grid, GridFittedBottom(bottom_height)) sea_ice = sea_ice_simulation(sea_ice_grid) - liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus # Set the ocean temperature and salinity set!(ocean.model, T=temperature_metadata[1], S=salinity_metadata[1]) From cf5dd4fd5bf685f17dfd77bcd1dd7796bc69c996 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 16:53:33 +0100 Subject: [PATCH 28/56] thermodynamics to ice_thermodynamics --- .../InterfaceComputations/component_interfaces.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 29713832f..f0002999e 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -233,7 +233,7 @@ end default_ai_temperature(sea_ice) = nothing function default_ai_temperature(sea_ice::SeaIceSimulation) - conductive_flux = sea_ice.model.thermodynamics.internal_heat_flux.parameters.flux + conductive_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux.parameters.flux return SkinTemperature(conductive_flux) end @@ -316,7 +316,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; sea_ice_properties = (reference_density = sea_ice_reference_density, heat_capacity = sea_ice_heat_capacity, freshwater_density = freshwater_density, - liquidus = sea_ice.model.thermodynamics.phase_transitions.liquidus, + liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus, temperature_units = sea_ice_temperature_units) net_momentum_fluxes = if sea_ice.model.dynamics isa Nothing @@ -372,7 +372,7 @@ function sea_ice_similarity_theory(sea_ice::SeaIceSimulation) # SkinTemperature. Also we need to pass the sea ice internal flux # The thickness and salinity need to be passed as well, # but the can be passed as state variables once we refactor the `StateValues` struct. - internal_flux = sea_ice.model.thermodynamics.internal_heat_flux + internal_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux interface_temperature_type = SkinTemperature(internal_flux) return SimilarityTheoryFluxes(; interface_temperature_type) end From 83d73e0d4a1ab9e62b25f87ecfc5db95e98cee55 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 17:07:56 +0100 Subject: [PATCH 29/56] try it out --- experiments/arctic_simulation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 8a48fd302..3a3675e6b 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -6,7 +6,7 @@ using Oceananigans.Units using Oceananigans.OrthogonalSphericalShellGrids using ClimaOcean.OceanSimulations using ClimaOcean.ECCO -using ClimaOcean.ECCO: all_ECCO_dates +using ClimaOcean.DataWrangling using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium using Printf @@ -163,7 +163,7 @@ run!(arctic) ##### version = ECCO4Monthly() -dates = all_ECCO_dates(version)[1:12] +dates = all_dates(version)[1:12] h_metadata = ECCOMetadata(:sea_ice_thickness; version, dates) ℵ_metadata = ECCOMetadata(:sea_ice_concentration; version, dates) From ddafd107bba8d3aeb03baf3ea5b2d59534f7ce95 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 17:10:25 +0100 Subject: [PATCH 30/56] woops --- .../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 f0002999e..83fb6859b 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -198,7 +198,7 @@ function atmosphere_sea_ice_interface(sea_ice::SeaIceSimulation, temperature_formulation, velocity_formulation) - interface_temperature = sea_ice.model.thermodynamics.top_surface_temperature + interface_temperature = sea_ice.model.ice_thermodynamics.top_surface_temperature return AtmosphereInterface(fluxes, ai_flux_formulation, interface_temperature, properties) end From 7395d24288faea688279f353e392c1c814de5d57 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 17:11:07 +0100 Subject: [PATCH 31/56] bring simulation up to speed --- experiments/arctic_simulation.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 3a3675e6b..5e4442dde 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -80,8 +80,8 @@ dynamics = SeaIceMomentumEquation(grid; sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) -set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), - ℵ=ECCOMetadata(:sea_ice_concentration)) +set!(sea_ice.model, h=Metadata(:sea_ice_thickness), + ℵ=Metadata(:sea_ice_concentration)) ##### ##### A Prescribed Atmosphere model @@ -162,11 +162,11 @@ run!(arctic) ##### Comparison to ECCO Climatology ##### -version = ECCO4Monthly() +dataset = ECCO4Monthly() dates = all_dates(version)[1:12] -h_metadata = ECCOMetadata(:sea_ice_thickness; version, dates) -ℵ_metadata = ECCOMetadata(:sea_ice_concentration; version, dates) +h_metadata = Metadata(:sea_ice_thickness; dataset, dates) +ℵ_metadata = Metadata(:sea_ice_concentration; dataset, dates) # Montly averaged ECCO data hE = ECCOFieldTimeSeries(h_metadata, grid; time_indices_in_memory=12) From f77da6041644ffd43cfe9531b3df9eb20b5f625a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 17:12:29 +0100 Subject: [PATCH 32/56] also here --- experiments/arctic_simulation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 5e4442dde..f11a0cbf8 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -49,8 +49,8 @@ ocean = ocean_simulation(grid; free_surface, closure) -set!(ocean.model, T=ECCOMetadata(:temperature), - S=ECCOMetadata(:salinity)) +set!(ocean.model, T=Metadata(:temperature), + S=Metadata(:salinity)) ##### ##### A Prognostic Sea-ice model From dd41cceee715517e1d6f0ea5f4dda473830a99c8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 11 Mar 2025 17:13:40 +0100 Subject: [PATCH 33/56] add a dataset --- experiments/arctic_simulation.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index f11a0cbf8..8af366219 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -49,8 +49,10 @@ ocean = ocean_simulation(grid; free_surface, closure) -set!(ocean.model, T=Metadata(:temperature), - S=Metadata(:salinity)) +dataset = ECCO4Monthly() + +set!(ocean.model, T=Metadata(:temperature; dataset), + S=Metadata(:salinity; dataset)) ##### ##### A Prognostic Sea-ice model @@ -80,8 +82,8 @@ dynamics = SeaIceMomentumEquation(grid; sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) -set!(sea_ice.model, h=Metadata(:sea_ice_thickness), - ℵ=Metadata(:sea_ice_concentration)) +set!(sea_ice.model, h=Metadata(:sea_ice_thickness; dataset), + ℵ=Metadata(:sea_ice_concentration; dataset)) ##### ##### A Prescribed Atmosphere model From addfdcb06a6c8be5d14dbb79f69ce6cfda034903 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 12 Mar 2025 09:51:39 +0100 Subject: [PATCH 34/56] correct similarity fluxes --- .../InterfaceComputations/component_interfaces.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 83fb6859b..284a541ef 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -266,9 +266,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)), - atmosphere_sea_ice_flux_formulation = CoefficientBasedFluxes(drag_coefficient=2e-3, - heat_transfer_coefficient=1e-4, - vapor_flux_coefficient=1e-4), + atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)), atmosphere_ocean_interface_temperature = BulkTemperature(), atmosphere_ocean_velocity_difference = RelativeVelocity(), atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean), From d4325dd7486e2a61bd80241c956ccbb586a7f31b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 12 Mar 2025 09:53:41 +0100 Subject: [PATCH 35/56] we were not capping! --- .../InterfaceComputations/interface_states.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 7488e349b..481ad766c 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -214,8 +214,8 @@ end h = Ψᵢ.h # Bottom temperature at the melting temperature - Tᵢ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) - Tᵢ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ) + Tₘ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) + Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ) Tₛ⁻ = Ψₛ.T #= @@ -227,7 +227,6 @@ end T★ = Tᵢ - Qₐ * h / k - # Under heating fluxes, cap surface temperature by melting temperature Tₘ = ℙᵢ.liquidus.freshwater_melting_temperature Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) @@ -241,6 +240,9 @@ end abs_ΔT = min(max_ΔT, abs(ΔT★)) Tₛ⁺ = Tₛ⁻ + abs_ΔT * sign(ΔT★) + # Under heating fluxes, cap surface temperature by melting temperature + Tₛ⁺ = min(Tₛ⁺, Tₘ) + return Tₛ⁺ end From 135296a313dbdce55736e0a5292ed64aaaf6bda3 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 12 Mar 2025 10:03:29 +0100 Subject: [PATCH 36/56] bugfix --- src/OceanSeaIceModels/InterfaceComputations/interface_states.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 481ad766c..0cabee73e 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -215,7 +215,7 @@ end # Bottom temperature at the melting temperature Tₘ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) - Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ) + Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) Tₛ⁻ = Ψₛ.T #= From 7c169f0a9d361cbeb2d47f51434c3f0611ce7638 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 12 Mar 2025 10:26:14 +0100 Subject: [PATCH 37/56] remove a mistake --- .../InterfaceComputations/interface_states.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 0cabee73e..108c9413e 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -214,8 +214,8 @@ end h = Ψᵢ.h # Bottom temperature at the melting temperature - Tₘ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) - Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) + Tᵢ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) + Tᵢ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ) Tₛ⁻ = Ψₛ.T #= @@ -227,9 +227,6 @@ end T★ = Tᵢ - Qₐ * h / k - Tₘ = ℙᵢ.liquidus.freshwater_melting_temperature - Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) - # Fix a NaN T★ = ifelse(isnan(T★), Tₛ⁻, T★) @@ -241,6 +238,8 @@ end Tₛ⁺ = Tₛ⁻ + abs_ΔT * sign(ΔT★) # Under heating fluxes, cap surface temperature by melting temperature + Tₘ = ℙᵢ.liquidus.freshwater_melting_temperature + Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) Tₛ⁺ = min(Tₛ⁺, Tₘ) return Tₛ⁺ From b757c7855820e8933d71417bcdf333e15606dcdf Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 12 Mar 2025 05:27:15 -0400 Subject: [PATCH 38/56] on gpus --- experiments/arctic_simulation.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 8af366219..6538434c9 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -12,6 +12,7 @@ using Printf using CUDA CUDA.device!(1) +arch = GPU() r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=2000) z_faces = MutableVerticalDiscretization(r_faces) @@ -20,13 +21,13 @@ Nx = 180 # longitudinal direction -> 250 points is about 1.5ᵒ resolution Ny = 180 # meridional direction -> same thing, 48 points is about 1.5ᵒ resolution Nz = length(r_faces) - 1 -grid = RotatedLatitudeLongitudeGrid(GPU(), size = (Nx, Ny, Nz), - latitude = (-45, 45), - longitude = (-45, 45), - z = r_faces, - north_pole = (180, 0), - halo = (5, 5, 4), - topology = (Bounded, Bounded, Bounded)) +grid = RotatedLatitudeLongitudeGrid(arch, size = (Nx, Ny, Nz), + latitude = (-45, 45), + longitude = (-45, 45), + z = r_faces, + north_pole = (180, 0), + halo = (5, 5, 4), + topology = (Bounded, Bounded, Bounded)) bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) @@ -89,7 +90,7 @@ set!(sea_ice.model, h=Metadata(:sea_ice_thickness; dataset), ##### A Prescribed Atmosphere model ##### -atmosphere = JRA55PrescribedAtmosphere(GPU(); backend=JRA55NetCDFBackend(40)) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40)) radiation = Radiation() ##### From d58a7c48a9419ddd19c892a19885d1f4428f4f73 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Thu, 13 Mar 2025 07:03:06 -0400 Subject: [PATCH 39/56] some kernel fusion --- .../atmosphere_sea_ice_fluxes.jl | 3 +- .../InterfaceComputations/interface_states.jl | 2 +- .../sea_ice_ocean_fluxes.jl | 80 +++++-------------- 3 files changed, 24 insertions(+), 61 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index b806460c5..c09e60f32 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -75,8 +75,7 @@ end i, j = @index(Global, NTuple) kᴺ = size(grid, 3) # index of the top ocean cell - time = Time(clock.time) - FT = eltype(grid) + FT = eltype(grid) @inbounds begin uₐ = atmosphere_state.u[i, j, 1] diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 108c9413e..8932a7e83 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -157,7 +157,7 @@ struct SkinTemperature{I, FT} max_ΔT :: FT end -SkinTemperature(internal_flux; max_ΔT=5) = SkinTemperature(internal_flux, max_ΔT) +SkinTemperature(internal_flux; max_ΔT=20) = SkinTemperature(internal_flux, max_ΔT) struct DiffusiveFlux{Z, K} δ :: Z # Boundary layer thickness, as a first guess we will use half the grid spacing diff --git a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index 72769e130..0ef926c9c 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -2,23 +2,19 @@ using Oceananigans.Operators: Δzᶜᶜᶜ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature function compute_sea_ice_ocean_fluxes!(coupled_model) - compute_sea_ice_ocean_salinity_flux!(coupled_model) - compute_sea_ice_ocean_latent_heat_flux!(coupled_model) - return nothing -end - -function compute_sea_ice_ocean_latent_heat_flux!(coupled_model) ocean = coupled_model.ocean sea_ice = coupled_model.sea_ice - Qᶠₒ = coupled_model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat - Qᵢₒ = coupled_model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat - + + sea_ice_ocean_fluxes = coupled_model.interfaces.sea_ice_ocean_interface.fluxes interface_properties = coupled_model.interfaces.sea_ice_ocean_interface.properties + Δt = ocean.Δt Tₒ = ocean.model.tracers.T Sₒ = ocean.model.tracers.S - Δt = ocean.Δt + Sᵢ = sea_ice.model.tracers.S ℵᵢ = sea_ice.model.ice_concentration + hᵢ = sea_ice.model.ice_thickness + h⁻ = coupled_model.interfaces.sea_ice_ocean_interface.previous_ice_thickness ocean_properties = coupled_model.interfaces.ocean_properties liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus @@ -28,15 +24,17 @@ function compute_sea_ice_ocean_latent_heat_flux!(coupled_model) # What about the latent heat removed from the ocean when ice forms? # Is it immediately removed from the ocean? Or is it stored in the ice? launch!(arch, grid, :xy, _compute_sea_ice_ocean_latent_heat_flux!, - Qᶠₒ, Qᵢₒ, grid, ℵᵢ, Tₒ, Sₒ, liquidus, ocean_properties, interface_properties, Δt) + sea_ice_ocean_fluxes, grid, hᵢ, h⁻, ℵᵢ, Sᵢ, Tₒ, Sₒ, liquidus, ocean_properties, interface_properties, Δt) return nothing end -@kernel function _compute_sea_ice_ocean_latent_heat_flux!(frazil_heat_flux, - interface_heat_flux, +@kernel function _compute_sea_ice_ocean_latent_heat_flux!(sea_ice_ocean_fluxes, grid, + ice_thickness, + previous_ice_thickness, ice_concentration, + ice_salinity, ocean_temperature, ocean_salinity, liquidus, @@ -47,10 +45,14 @@ end i, j = @index(Global, NTuple) Nz = size(grid, 3) - Qᶠₒ = frazil_heat_flux - Qᵢₒ = interface_heat_flux + Qᶠₒ = sea_ice_ocean_fluxes.frazil_heat + Qᵢₒ = sea_ice_ocean_fluxes.interface_heat + Jˢ = sea_ice_ocean_fluxes.salt Tₒ = ocean_temperature Sₒ = ocean_salinity + Sᵢ = ice_salinity + hᵢ = ice_thickness + h⁻ = previous_ice_thickness ρₒ = ocean_properties.reference_density cₒ = ocean_properties.heat_capacity uₘ★ = interface_properties.characteristic_melting_speed @@ -102,54 +104,16 @@ end # A positive value δQ_melting > 0 corresponds to ocean cooling; ie # is fluxing upwards, into the ice. This occurs when applying the # ice bath equilibrium condition to cool down a warm ocean (δEₒ < 0). - δQ_melting = - δE_ice_bath * uₘ★ + Δz = Δzᶜᶜᶜ(i, j, Nz, grid) + δQ_melting = - δE_ice_bath * uₘ★ # Store column-integrated ice-ocean heat flux @inbounds Qᶠₒ[i, j, 1] = δQ_frazil @inbounds Qᵢₒ[i, j, 1] = δQ_melting * ℵ # Melting depends on concentration -end - -function compute_sea_ice_ocean_salinity_flux!(coupled_model) - # Compute salinity increment due to changes in ice thickness - - sea_ice = coupled_model.sea_ice - ocean = coupled_model.ocean - grid = sea_ice.model.grid - arch = architecture(grid) - Sₒ = ocean.model.tracers.S - Sᵢ = sea_ice.model.tracers.S - Δt = ocean.Δt - hⁿ = sea_ice.model.ice_thickness - h⁻ = coupled_model.interfaces.sea_ice_ocean_interface.previous_ice_thickness - - interface_fluxes = coupled_model.interfaces.sea_ice_ocean_interface.fluxes - - launch!(arch, grid, :xy, _compute_sea_ice_ocean_salinity_flux!, - interface_fluxes.salt, grid, hⁿ, h⁻, Sᵢ, Sₒ, Δt) - - return nothing -end - -@kernel function _compute_sea_ice_ocean_salinity_flux!(salt_flux, - grid, - ice_thickness, - previous_ice_thickness, - ice_salinity, - ocean_salinity, - Δt) - i, j = @index(Global, NTuple) - - Nz = size(grid, 3) - - hⁿ = ice_thickness - h⁻ = previous_ice_thickness - Sᵢ = ice_salinity - Sₒ = ocean_salinity - Jˢ = salt_flux @inbounds begin # Change in thickness - Δh = hⁿ[i, j, 1] - h⁻[i, j, 1] + Δh = hᵢ[i, j, 1] - h⁻[i, j, 1] # Update surface salinity flux. # Note: the Δt below is the ocean time-step, eg. @@ -157,6 +121,6 @@ end Jˢ[i, j, 1] = Δh / Δt * (Sᵢ[i, j, 1] - Sₒ[i, j, Nz]) # Update previous ice thickness - h⁻[i, j, 1] = hⁿ[i, j, 1] + h⁻[i, j, 1] = hᵢ[i, j, 1] end -end +end \ No newline at end of file From ceda2174dd6be689abab0ca4fc013f2507dbabc0 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Thu, 13 Mar 2025 07:06:24 -0400 Subject: [PATCH 40/56] back to how it was --- src/OceanSeaIceModels/InterfaceComputations/interface_states.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 8932a7e83..108c9413e 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -157,7 +157,7 @@ struct SkinTemperature{I, FT} max_ΔT :: FT end -SkinTemperature(internal_flux; max_ΔT=20) = SkinTemperature(internal_flux, max_ΔT) +SkinTemperature(internal_flux; max_ΔT=5) = SkinTemperature(internal_flux, max_ΔT) struct DiffusiveFlux{Z, K} δ :: Z # Boundary layer thickness, as a first guess we will use half the grid spacing From 020a47a5a4eeb6f76be4238ceda793d8caafec3f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 13:48:27 +0100 Subject: [PATCH 41/56] add this --- experiments/single_column_sea_ice.jl | 184 +++++++++++++++++++++++++++ 1 file changed, 184 insertions(+) create mode 100644 experiments/single_column_sea_ice.jl diff --git a/experiments/single_column_sea_ice.jl b/experiments/single_column_sea_ice.jl new file mode 100644 index 000000000..dea799a13 --- /dev/null +++ b/experiments/single_column_sea_ice.jl @@ -0,0 +1,184 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf + +grid = RectilinearGrid(arch, size = (), + latitude = 87.0, + longitude = 33.0, + topology = (Flat, Flat, Flat)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) + +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) + +##### +##### A Propgnostic Ocean model +##### + +# A very diffusive ocean +momentum_advection = WENOVectorInvariant(order=3) +tracer_advection = WENO(order=3) + +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) +closure = ClimaOcean.OceanSimulations.default_ocean_closure() + +ocean = ocean_simulation(grid; + momentum_advection, + tracer_advection, + free_surface, + closure) + +dataset = ECCO4Monthly() + +set!(ocean.model, T=Metadata(:temperature; dataset), + S=Metadata(:salinity; dataset)) + +##### +##### A Prognostic Sea-ice model +##### + +using ClimaSeaIce.SeaIceMomentumEquations +using ClimaSeaIce.Rheologies + +# Remember to pass the SSS as a bottom bc to the sea ice! +SSS = view(ocean.model.tracers.S, :, :, grid.Nz) +bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS) + +SSU = view(ocean.model.velocities.u, :, :, grid.Nz) +SSV = view(ocean.model.velocities.u, :, :, grid.Nz) + +τo = SemiImplicitStress(uₑ=SSU, vₑ=SSV) +τua = Field{Face, Center, Nothing}(grid) +τva = Field{Center, Face, Nothing}(grid) + +dynamics = SeaIceMomentumEquation(grid; + coriolis = ocean.model.coriolis, + top_momentum_stress = (u=τua, v=τva), + bottom_momentum_stress = τo, + ocean_velocities = (u=SSU, v=SSV), + rheology = ElastoViscoPlasticRheology(), + solver = SplitExplicitSolver(120)) + +sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) + +set!(sea_ice.model, h=Metadata(:sea_ice_thickness; dataset), + ℵ=Metadata(:sea_ice_concentration; dataset)) + +##### +##### A Prescribed Atmosphere model +##### + +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40)) +radiation = Radiation() + +##### +##### Arctic coupled model +##### + +arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) + +# Sea-ice variables +h = sea_ice.model.ice_thickness +ℵ = sea_ice.model.ice_concentration +u = sea_ice.model.velocities.u +v = sea_ice.model.velocities.v + +# Fluxes +Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature +Qˡ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat +Qˢ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat +Qⁱ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat +Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat +Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat +Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat +τx = arctic.model.interfaces.net_fluxes.sea_ice_top.u +τy = arctic.model.interfaces.net_fluxes.sea_ice_top.v + +# Output writers +arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, u, v, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, τx, τy), + filename = "sea_ice_quantities.jld2", + schedule = IterationInterval(12), + overwrite_existing=true) + +arctic.output_writers[:averages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, u, v, τx, τy), + filename = "averaged_sea_ice_quantities.jld2", + schedule = AveragedTimeInterval(1days), + overwrite_existing=true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + hmean = mean(sea_ice.model.ice_thickness) + ℵmean = mean(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg3 = @sprintf("mean(h): %.2e m, mean(ℵ): %.2e ", hmean, ℵmean) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 * msg4 * msg5 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(arctic, progress, IterationInterval(10)) + +run!(arctic) + +##### +##### Comparison to ECCO Climatology +##### + +dataset = ECCO4Monthly() +dates = all_dates(version)[1:12] + +h_metadata = Metadata(:sea_ice_thickness; dataset, dates) +ℵ_metadata = Metadata(:sea_ice_concentration; dataset, dates) + +# Montly averaged ECCO data +hE = ECCOFieldTimeSeries(h_metadata, grid; time_indices_in_memory=12) +ℵE = ECCOFieldTimeSeries(ℵ_metadata, grid; time_indices_in_memory=12) + +# Daily averaged Model output +h = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "h") +ℵ = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "ℵ") + +# Montly average the model output +hm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) +ℵm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) + +for (i, time) in enumerate(hm.times) + counter = 0 + for (j, t) in enumerate(h.times) + if t ≤ time + hm[i] .+= h[j] + ℵm[i] .+= ℵ[j] + counter += 1 + end + end + @show counter + hm[i] ./= counter + ℵm[i] ./= counter +end From 16655a26183cf32258ea342f4c4558e66326a7fd Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:04:09 +0100 Subject: [PATCH 42/56] Revert "Revert "add a prescribed ocean"" This reverts commit 71d02b840dedf8378d43429af23a709d7d793e8b. --- .../flux_climatology/flux_climatology.jl | 55 ------------- src/ClimaOcean.jl | 10 +++ src/OceanSeaIceModels/OceanSeaIceModels.jl | 1 - src/OceanSeaIceModels/ocean_sea_ice_model.jl | 19 ----- src/OceanSimulations/OceanSimulations.jl | 12 ++- .../prescribed_ocean_model.jl | 81 +++++++++++++++++++ src/SeaIceSimulations.jl | 5 ++ test/test_ocean_sea_ice_model.jl | 19 +++++ 8 files changed, 126 insertions(+), 76 deletions(-) create mode 100644 src/OceanSimulations/prescribed_ocean_model.jl diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 94a1a951e..f0f766663 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -115,32 +115,6 @@ end ##### A prescribed ocean... ##### -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} - architecture :: A - grid :: G - clock :: Clock{C} - velocities :: U - tracers :: T - timeseries :: F -end - -function PrescribedOcean(timeseries; - grid, - clock=Clock{Float64}(time = 0)) - - τx = Field{Face, Center, Nothing}(grid) - τy = Field{Center, Face, Nothing}(grid) - Jᵀ = Field{Center, Center, Nothing}(grid) - Jˢ = Field{Center, Center, Nothing}(grid) - - u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) - v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) - T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) - S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - - PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) -end - # ...with prescribed velocity and tracer fields version = ECCO4Monthly() dates = all_ECCO_dates(version)[1:24] @@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch) ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days) -##### -##### Need to extend a couple of methods -##### - -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) - tick!(model.clock, Δt) - time = Time(model.clock.time) - - possible_fts = merge(model.velocities, model.tracers) - time_series_tuple = extract_field_time_series(possible_fts) - - for fts in time_series_tuple - update_field_time_series!(fts, time) - end - - parent(model.velocities.u) .= parent(model.timeseries.u[time]) - parent(model.velocities.v) .= parent(model.timeseries.v[time]) - parent(model.tracers.T) .= parent(model.timeseries.T[time]) - parent(model.tracers.S) .= parent(model.timeseries.S[time]) - - return nothing -end - -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing - -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 - ##### ##### A prescribed atmosphere... ##### diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 2a90d73e8..96a543673 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -41,6 +41,16 @@ using DataDeps using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries using Oceananigans.Grids: node +# Does not really matter if there is no model +reference_density(::Nothing) = 0 +heat_capacity(::Nothing) = 0 + +reference_density(unsupported) = + throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))")) + +heat_capacity(unsupported) = + throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))")) + const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries, GPUAdaptedFieldTimeSeries} diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 8041f497d..4dee01e38 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -22,7 +22,6 @@ using ClimaSeaIce: SeaIceModel using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using ClimaOcean: stateindex - using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index d2344dabb..9c97a6cb5 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -73,25 +73,6 @@ function initialize!(model::OSIM) return nothing end -reference_density(unsupported) = - throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))")) - -heat_capacity(unsupported) = - throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))")) - -reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation) -reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) -reference_density(eos::TEOS10EquationOfState) = eos.reference_density -reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density - -heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation) -heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) -heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity - -# Does not really matter if there is no model -reference_density(::Nothing) = 0 -heat_capacity(::Nothing) = 0 - function heat_capacity(::TEOS10EquationOfState{FT}) where FT cₚ⁰ = SeawaterPolynomials.TEOS10.teos10_reference_heat_capacity return convert(FT, cₚ⁰) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 32423c690..49d48aa6e 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation +export ocean_simulation, PrescribedOceanModel using Oceananigans using Oceananigans.Units @@ -21,6 +21,15 @@ using Oceananigans.BuoyancyFormulations: g_Earth using Oceananigans.Coriolis: Ω_Earth using Oceananigans.Operators +import ClimaOcean: reference_density, heat_capacity + +reference_density(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = reference_density(ocean.model.buoyancy.formulation) +reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state) +reference_density(eos::TEOS10EquationOfState) = eos.reference_density + +heat_capacity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = heat_capacity(ocean.model.buoyancy.formulation) +heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state) + struct Default{V} value :: V end @@ -41,5 +50,6 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("ocean_simulation.jl") +include("prescribed_ocean_model.jl") end # module diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl new file mode 100644 index 000000000..a5309e4a8 --- /dev/null +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -0,0 +1,81 @@ +using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! + +import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! +import Oceananigans.Models: timestepper, update_model_field_time_series! + +import ClimaOcean: reference_density, heat_capacity +import Oceananigans.Architectures: on_architecture + +##### +##### A prescribed ocean... +##### + +struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} + architecture :: A + grid :: G + clock :: Clock{C} + velocities :: U + tracers :: T + timeseries :: F +end + +function PrescribedOcean(timeseries; + grid, + clock=Clock{Float64}(time = 0)) + + # Make sure all elements of the timeseries are on the same grid + for (k, v) in timeseries + if !isa(v, FieldTimeSeries) + throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) + end + if v.grid != grid + throw(ArgumentError("All variables in the timeseries reside on the provided grid")) + end + end + + τx = Field{Face, Center, Nothing}(grid) + τy = Field{Center, Face, Nothing}(grid) + Jᵀ = Field{Center, Center, Nothing}(grid) + Jˢ = Field{Center, Center, Nothing}(grid) + + u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx))) + v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy))) + T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) + S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) + + PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) +end + +##### +##### Need to extend a couple of methods +##### + +function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) + tick!(model.clock, Δt) + time = Time(model.clock.time) + + possible_fts = merge(model.velocities, model.tracers) + time_series_tuple = extract_field_time_series(possible_fts) + + for fts in time_series_tuple + update_field_time_series!(fts, time) + end + + update_u_velocity = haskey(model.timeseries, :u) + update_v_velocity = haskey(model.timeseries, :v) + update_temperature = haskey(model.timeseries, :T) + update_salinity = haskey(model.timeseries, :S) + + update_u_velocity && iterpolate!(model.velocities.u, model.timeseries.u[time]) + update_v_velocity && iterpolate!(model.velocities.v, model.timeseries.v[time]) + update_temperature && iterpolate!(model.tracers.T, model.timeseries.T[time]) + update_salinity && iterpolate!(model.tracers.S, model.timeseries.S[time]) + + return nothing +end + +update_state!(::PrescribedOcean) = nothing +timestepper(::PrescribedOcean) = nothing + +reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 2b7dc17b8..510680651 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -20,6 +20,11 @@ using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium using ClimaOcean.OceanSimulations: Default +import ClimaOcean: reference_density, heat_capacity + +reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density +heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity + function sea_ice_simulation(grid; Δt = 5minutes, ice_salinity = 0, # psu diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl index 55b6ac3f5..16d55fc9d 100644 --- a/test/test_ocean_sea_ice_model.jl +++ b/test/test_ocean_sea_ice_model.jl @@ -100,5 +100,24 @@ using ClimaSeaIce.Rheologies true end end + + @testset "Test a prescribed ocean model" begin + grid = LatitudeLongitudeGrid(size = (10, 10, 10), latitude = (20, 30), longitude = (40, 50), z = (-100, 0)) + T = ECCOFieldTimeSeries(:temperature, grid; time_indices_in_memory = 2) + S = ECCOFieldTimeSeries(:salinity, grid; time_indices_in_memory = 2) + ocean_model = PrescribedOcean((; T, S); grid) + + time_step!(ocean_model, T.times[2]) + + @test ocean_model.clock.time == T.times[2] + @test ocean_model.tracers.T.data == T[2].data + @test ocean_model.tracers.S.data == S[2].data + + time_step!(ocean_model, 10days) + + @test ocean_model.clock.time == 10days + @test T.backend.start == 2 + @test ocean_model.tracers.T.data != T[3].data + end end From 1a00396ddb397a16107817590427bcc3ad0ad5c2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:09:13 +0100 Subject: [PATCH 43/56] better like this --- .../prescribed_ocean_model.jl | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index a5309e4a8..e2bf596fb 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -19,17 +19,18 @@ struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} timeseries :: F end -function PrescribedOcean(timeseries; +function PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) # Make sure all elements of the timeseries are on the same grid - for (k, v) in timeseries - if !isa(v, FieldTimeSeries) - throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) - end - if v.grid != grid - throw(ArgumentError("All variables in the timeseries reside on the provided grid")) + # If we decide to pass a timeseries + if !isempty(timeseries) + for (k, v) in timeseries + isa(v, FieldTimeSeries) || + throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) + v.grid == grid || + throw(ArgumentError("All variables in the timeseries reside on the provided grid")) end end @@ -66,10 +67,10 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) update_temperature = haskey(model.timeseries, :T) update_salinity = haskey(model.timeseries, :S) - update_u_velocity && iterpolate!(model.velocities.u, model.timeseries.u[time]) - update_v_velocity && iterpolate!(model.velocities.v, model.timeseries.v[time]) - update_temperature && iterpolate!(model.tracers.T, model.timeseries.T[time]) - update_salinity && iterpolate!(model.tracers.S, model.timeseries.S[time]) + update_u_velocity && parent(model.velocities.u) .= parent(model.timeseries.u[time]) + update_v_velocity && parent(model.velocities.v) .= parent(model.timeseries.v[time]) + update_temperature && parent(model.tracers.T) .= parent(model.timeseries.T[time]) + update_salinity && parent(model.tracers.S) .= parent(model.timeseries.S[time]) return nothing end From 7e5babd0c494a741babc0b49170b45bd5a27b802 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:11:37 +0100 Subject: [PATCH 44/56] add a docstring --- src/OceanSimulations/prescribed_ocean_model.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index e2bf596fb..5d0c42792 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -10,8 +10,8 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} - architecture :: A +struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} + architecture :: Arch grid :: G clock :: Clock{C} velocities :: U @@ -19,6 +19,18 @@ struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} timeseries :: F end +""" + PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + +Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` +on a `grid` with a `clock`. + +Arguments +========= +- `timeseries::NamedTuple`: A named tuple containing time series data for various fields. The named tuple can be empty or + include any combination of the following fields: `u`, `v`, `T`, `S`. + All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. +""" function PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) @@ -44,7 +56,7 @@ function PrescribedOcean(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) end ##### From 83b1a8ab8df52551d8c5f0821517b84d9b2e3c8c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 14:12:24 +0100 Subject: [PATCH 45/56] better comment --- src/OceanSimulations/prescribed_ocean_model.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 5d0c42792..15a73ca33 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -27,9 +27,10 @@ on a `grid` with a `clock`. Arguments ========= -- `timeseries::NamedTuple`: A named tuple containing time series data for various fields. The named tuple can be empty or - include any combination of the following fields: `u`, `v`, `T`, `S`. - All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. +- `timeseries`: A `NamedTuple` containing time series data for various fields. The named tuple can be empty + (meaning that the ocean does not evolve in time) or include any combination of the + following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` + and reside on the provided `grid`. """ function PrescribedOcean(timeseries=NamedTuple(); grid, From 06db8a3ba9d8916419119bd9f187fe4cf23a678a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:26:47 +0100 Subject: [PATCH 46/56] arctic single column model --- experiments/single_column_sea_ice.jl | 145 +++++------------- src/DataWrangling/ECCO/ECCO_restoring.jl | 5 +- .../InterfaceComputations/radiation.jl | 2 +- src/OceanSeaIceModels/OceanSeaIceModels.jl | 4 +- .../prescribed_ocean_model.jl | 51 +++--- src/SeaIceSimulations.jl | 4 +- 6 files changed, 77 insertions(+), 134 deletions(-) diff --git a/experiments/single_column_sea_ice.jl b/experiments/single_column_sea_ice.jl index dea799a13..092197f18 100644 --- a/experiments/single_column_sea_ice.jl +++ b/experiments/single_column_sea_ice.jl @@ -3,94 +3,71 @@ using ClimaSeaIce using Oceananigans using Oceananigans.Grids using Oceananigans.Units -using Oceananigans.OrthogonalSphericalShellGrids using ClimaOcean.OceanSimulations using ClimaOcean.ECCO using ClimaOcean.DataWrangling using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium using Printf -grid = RectilinearGrid(arch, size = (), +# A single column grid with one point in the Arctic ocean +grid = LatitudeLongitudeGrid(size = 1, latitude = 87.0, longitude = 33.0, - topology = (Flat, Flat, Flat)) - -bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) - -grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) + z = (-10, 0), + topology = (Flat, Flat, Bounded)) ##### -##### A Propgnostic Ocean model +##### A Prescribed Ocean model ##### -# A very diffusive ocean -momentum_advection = WENOVectorInvariant(order=3) -tracer_advection = WENO(order=3) +dataset = ECCO4Monthly() +dates = all_dates(dataset)[1:24] -free_surface = SplitExplicitFreeSurface(grid; cfl=0.7) -closure = ClimaOcean.OceanSimulations.default_ocean_closure() +T = ECCOFieldTimeSeries(:temperature, grid; dates, inpainting=nothing, time_indices_in_memory=24) +S = ECCOFieldTimeSeries(:salinity, grid; dates, inpainting=nothing, time_indices_in_memory=24) -ocean = ocean_simulation(grid; - momentum_advection, - tracer_advection, - free_surface, - closure) +ocean_model = PrescribedOceanModel((; T, S); grid) +ocean = Simulation(ocean_model, Δt=30minutes, stop_time=730days) -dataset = ECCO4Monthly() +##### +##### A Prescribed Atmosphere model +##### -set!(ocean.model, T=Metadata(:temperature; dataset), - S=Metadata(:salinity; dataset)) +atmosphere = JRA55PrescribedAtmosphere(latitude = 87.0, + longitude = 33.0, + backend = InMemory()) ##### -##### A Prognostic Sea-ice model +##### A Prognostic Sea-ice model with no dynamics ##### -using ClimaSeaIce.SeaIceMomentumEquations -using ClimaSeaIce.Rheologies - # Remember to pass the SSS as a bottom bc to the sea ice! -SSS = view(ocean.model.tracers.S, :, :, grid.Nz) +SSS = view(ocean.model.tracers.S, :, :, 1) bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS) -SSU = view(ocean.model.velocities.u, :, :, grid.Nz) -SSV = view(ocean.model.velocities.u, :, :, grid.Nz) - -τo = SemiImplicitStress(uₑ=SSU, vₑ=SSV) -τua = Field{Face, Center, Nothing}(grid) -τva = Field{Center, Face, Nothing}(grid) - -dynamics = SeaIceMomentumEquation(grid; - coriolis = ocean.model.coriolis, - top_momentum_stress = (u=τua, v=τva), - bottom_momentum_stress = τo, - ocean_velocities = (u=SSU, v=SSV), - rheology = ElastoViscoPlasticRheology(), - solver = SplitExplicitSolver(120)) - -sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition) -set!(sea_ice.model, h=Metadata(:sea_ice_thickness; dataset), - ℵ=Metadata(:sea_ice_concentration; dataset)) +set!(sea_ice.model, h=Metadata(:sea_ice_thickness; dataset, dates=first(dates)), + ℵ=Metadata(:sea_ice_concentration; dataset, dates=first(dates))) -##### -##### A Prescribed Atmosphere model +##### +##### Radiation of the Ocean and the Sea ice ##### -atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40)) -radiation = Radiation() +# to mimick snow, as in Semtner (1975) we use a sea ice albedo that is 0.64 for top temperatures +# above -0.1 ᵒC and 0.75 for top temperatures below -0.1 ᵒC. We also disable radiation for the ocean. +radiation = Radiation(sea_ice_albedo=0.7, ocean_albedo=1, ocean_emissivity=0) ##### ##### Arctic coupled model ##### arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) +arctic = Simulation(arctic, Δt=30minutes, stop_time=730days) # Sea-ice variables h = sea_ice.model.ice_thickness ℵ = sea_ice.model.ice_concentration -u = sea_ice.model.velocities.u -v = sea_ice.model.velocities.v # Fluxes Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature @@ -100,16 +77,14 @@ Qⁱ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat -τx = arctic.model.interfaces.net_fluxes.sea_ice_top.u -τy = arctic.model.interfaces.net_fluxes.sea_ice_top.v # Output writers -arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, u, v, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, τx, τy), +arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "sea_ice_quantities.jld2", schedule = IterationInterval(12), overwrite_existing=true) -arctic.output_writers[:averages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, u, v, τx, τy), +arctic.output_writers[:averages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), filename = "averaged_sea_ice_quantities.jld2", schedule = AveragedTimeInterval(1days), overwrite_existing=true) @@ -120,65 +95,23 @@ using Statistics function progress(sim) sea_ice = sim.model.sea_ice - hmax = maximum(sea_ice.model.ice_thickness) - ℵmax = maximum(sea_ice.model.ice_concentration) - hmean = mean(sea_ice.model.ice_thickness) - ℵmean = mean(sea_ice.model.ice_concentration) - Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) - Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + h = first(sea_ice.model.ice_thickness) + ℵ = first(sea_ice.model.ice_concentration) + T = first(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) step_time = 1e-9 * (time_ns() - wall_time[]) msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) - msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) - msg3 = @sprintf("mean(h): %.2e m, mean(ℵ): %.2e ", hmean, ℵmean) - msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) - msg5 = @sprintf("wall time: %s \n", prettytime(step_time)) + msg2 = @sprintf("h: %.2e m, ℵ: %.2e , T: %.2e ᴼC", h, ℵ, T) + msg3 = @sprintf("wall time: %s \n", prettytime(step_time)) - @info msg1 * msg2 * msg3 * msg4 * msg5 + @info msg1 * msg2 * msg3 - wall_time[] = time_ns() - - return nothing + wall_time[] = time_ns() + return nothing end # And add it as a callback to the simulation. add_callback!(arctic, progress, IterationInterval(10)) -run!(arctic) - -##### -##### Comparison to ECCO Climatology -##### - -dataset = ECCO4Monthly() -dates = all_dates(version)[1:12] - -h_metadata = Metadata(:sea_ice_thickness; dataset, dates) -ℵ_metadata = Metadata(:sea_ice_concentration; dataset, dates) - -# Montly averaged ECCO data -hE = ECCOFieldTimeSeries(h_metadata, grid; time_indices_in_memory=12) -ℵE = ECCOFieldTimeSeries(ℵ_metadata, grid; time_indices_in_memory=12) - -# Daily averaged Model output -h = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "h") -ℵ = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "ℵ") - -# Montly average the model output -hm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) -ℵm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory()) - -for (i, time) in enumerate(hm.times) - counter = 0 - for (j, t) in enumerate(h.times) - if t ≤ time - hm[i] .+= h[j] - ℵm[i] .+= ℵ[j] - counter += 1 - end - end - @show counter - hm[i] ./= counter - ℵm[i] ./= counter -end +run!(arctic) \ No newline at end of file diff --git a/src/DataWrangling/ECCO/ECCO_restoring.jl b/src/DataWrangling/ECCO/ECCO_restoring.jl index 23483d032..f76d195d1 100644 --- a/src/DataWrangling/ECCO/ECCO_restoring.jl +++ b/src/DataWrangling/ECCO/ECCO_restoring.jl @@ -138,15 +138,14 @@ function ECCOFieldTimeSeries(metadata::ECCOMetadata, grid::AbstractGrid; return fts end -function ECCOFieldTimeSeries(variable_name::Symbol; +function ECCOFieldTimeSeries(variable_name::Symbol, arch_or_grid=CPU(); dataset = ECCO4Monthly(), - architecture = CPU(), dates = all_dates(dataset, variable_name), dir = download_ECCO_cache, kw...) metadata = Metadata(variable_name, dates, dataset, dir) - return ECCOFieldTimeSeries(metadata, architecture; kw...) + return ECCOFieldTimeSeries(metadata, arch_or_grid; kw...) end # Variable names for restorable data diff --git a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl index faf17481c..551b8c4ed 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl @@ -67,7 +67,7 @@ function Base.show(io::IO, r::Radiation) print(io, "├── stefan_boltzmann_constant: ", prettysummary(σ), '\n') print(io, "├── emission: ", summary(r.emission), '\n') print(io, "│ ├── ocean: ", prettysummary(r.emission.ocean), '\n') - print(io, "│ └── sea_ice: ", prettysummary(r.emission.ocean), '\n') + print(io, "│ └── sea_ice: ", prettysummary(r.emission.sea_ice), '\n') print(io, "└── reflection: ", summary(r.reflection), '\n') print(io, " ├── ocean: ", prettysummary(r.reflection.ocean), '\n') print(io, " └── sea_ice: ", prettysummary(r.reflection.sea_ice)) diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 4dee01e38..8316a157f 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -25,10 +25,10 @@ using ClimaOcean: stateindex using KernelAbstractions: @kernel, @index using KernelAbstractions.Extras.LoopInfo: @unroll +import ClimaOcean: reference_density, heat_capacity + function downwelling_radiation end function freshwater_flux end -function reference_density end -function heat_capacity end const default_gravitational_acceleration = 9.80665 const default_freshwater_density = 1000 diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean_model.jl index 15a73ca33..e1ff82296 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -1,3 +1,5 @@ +using Oceananigans.Models: AbstractModel +using Oceananigans.Fields: ZeroField using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series! import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick! @@ -10,7 +12,7 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} +struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} architecture :: Arch grid :: G clock :: Clock{C} @@ -20,7 +22,7 @@ struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} end """ - PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` on a `grid` with a `clock`. @@ -32,17 +34,18 @@ Arguments following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. """ -function PrescribedOcean(timeseries=NamedTuple(); +function PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) # Make sure all elements of the timeseries are on the same grid # If we decide to pass a timeseries if !isempty(timeseries) - for (k, v) in timeseries - isa(v, FieldTimeSeries) || + for k in keys(timeseries) + f = timeseries[k] + isa(f, FieldTimeSeries) || throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`")) - v.grid == grid || + f.grid == grid || throw(ArgumentError("All variables in the timeseries reside on the provided grid")) end end @@ -57,14 +60,14 @@ function PrescribedOcean(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) end ##### ##### Need to extend a couple of methods ##### -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) +function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) tick!(model.clock, Δt) time = Time(model.clock.time) @@ -75,21 +78,29 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) update_field_time_series!(fts, time) end - update_u_velocity = haskey(model.timeseries, :u) - update_v_velocity = haskey(model.timeseries, :v) - update_temperature = haskey(model.timeseries, :T) - update_salinity = haskey(model.timeseries, :S) + # Time stepping the model! - update_u_velocity && parent(model.velocities.u) .= parent(model.timeseries.u[time]) - update_v_velocity && parent(model.velocities.v) .= parent(model.timeseries.v[time]) - update_temperature && parent(model.tracers.T) .= parent(model.timeseries.T[time]) - update_salinity && parent(model.tracers.S) .= parent(model.timeseries.S[time]) + if haskey(model.timeseries, :u) + arent(model.velocities.u) .= parent(model.timeseries.u[time]) + end + + if haskey(model.timeseries, :v) + parent(model.velocities.v) .= parent(model.timeseries.v[time]) + end + + if haskey(model.timeseries, :T) + parent(model.tracers.T) .= parent(model.timeseries.T[time]) + end + + if haskey(model.timeseries, :S) + parent(model.tracers.S) .= parent(model.timeseries.S[time]) + end return nothing end -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing +update_state!(::PrescribedOceanModel) = nothing +timestepper(::PrescribedOceanModel) = nothing -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 +reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6 diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 510680651..59576aa0c 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -22,8 +22,8 @@ using ClimaOcean.OceanSimulations: Default import ClimaOcean: reference_density, heat_capacity -reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density -heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity +reference_density(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density +heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity function sea_ice_simulation(grid; Δt = 5minutes, From d9131bbcc198ac94fbd0185cf59e4cfc589843ba Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:42:56 +0100 Subject: [PATCH 47/56] start now --- experiments/single_column_sea_ice.jl | 9 +++------ .../InterfaceComputations/atmosphere_ocean_fluxes.jl | 6 ++---- src/OceanSeaIceModels/ocean_sea_ice_model.jl | 8 ++------ 3 files changed, 7 insertions(+), 16 deletions(-) diff --git a/experiments/single_column_sea_ice.jl b/experiments/single_column_sea_ice.jl index 092197f18..63cc91d1d 100644 --- a/experiments/single_column_sea_ice.jl +++ b/experiments/single_column_sea_ice.jl @@ -84,11 +84,6 @@ arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Q schedule = IterationInterval(12), overwrite_existing=true) -arctic.output_writers[:averages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ), - filename = "averaged_sea_ice_quantities.jld2", - schedule = AveragedTimeInterval(1days), - overwrite_existing=true) - wall_time = Ref(time_ns()) using Statistics @@ -114,4 +109,6 @@ end # And add it as a callback to the simulation. add_callback!(arctic, progress, IterationInterval(10)) -run!(arctic) \ No newline at end of file +run!(arctic) + +h = FieldTimeSeries("sea_ice_quantities.jld2", "h") \ No newline at end of file diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl index 97627ceab..4a1b2c830 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl @@ -1,5 +1,6 @@ using Oceananigans.Operators: intrinsic_vector using Oceananigans.Grids: inactive_node +using ClimaOcean.OceanSimulations using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: thermodynamics_parameters, surface_layer_height, boundary_layer_height @@ -11,10 +12,7 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) arch = architecture(grid) clock = coupled_model.clock - ocean_state = (u = ocean.model.velocities.u, - v = ocean.model.velocities.v, - T = ocean.model.tracers.T, - S = ocean.model.tracers.S) + ocean_state = ocean_state(ocean) atmosphere_fields = coupled_model.interfaces.exchanger.exchange_atmosphere_state diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 9c97a6cb5..418ec41fe 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -141,19 +141,15 @@ function default_nan_checker(model::OceanSeaIceModel) return nan_checker end -@kernel function _above_freezing_ocean_temperature!(T, grid, S, ℵ, liquidus) +@kernel function _above_freezing_ocean_temperature!(T, grid, S, liquidus) i, j = @index(Global, NTuple) Nz = size(grid, 3) @inbounds begin - for k in 1:Nz-1 + for k in 1:Nz Tm = melting_temperature(liquidus, S[i, j, k]) T[i, j, k] = max(T[i, j, k], Tm) end - - ℵi = ℵ[i, j, 1] - Tm = melting_temperature(liquidus, S[i, j, Nz]) - T[i, j, Nz] = ifelse(ℵi > 0, Tm, T[i, j, Nz]) end end From 2550442921441a861d8c6d0c89c530f43e102a3a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:53:04 +0100 Subject: [PATCH 48/56] start with this... --- .../component_interfaces.jl | 4 ++-- src/OceanSimulations/OceanSimulations.jl | 12 ++++++++++-- ...ibed_ocean_model.jl => prescribed_ocean.jl} | 18 +++++++++--------- 3 files changed, 21 insertions(+), 13 deletions(-) rename src/OceanSimulations/{prescribed_ocean_model.jl => prescribed_ocean.jl} (83%) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 284a541ef..9b0f49fd4 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -245,7 +245,7 @@ function default_ao_specific_humidity(ocean) end """ - ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; + ComponentInterfaces(atmosphere, ocean, sea_ice; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(), @@ -262,7 +262,7 @@ end sea_ice_reference_density = reference_density(sea_ice), sea_ice_heat_capacity = heat_capacity(sea_ice)) """ -function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; +function ComponentInterfaces(atmosphere, ocean, sea_ice; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)), diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 49d48aa6e..353160a04 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation, PrescribedOceanModel +export ocean_simulation, PrescribedOcean, ocean_state using Oceananigans using Oceananigans.Units @@ -50,6 +50,14 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("ocean_simulation.jl") -include("prescribed_ocean_model.jl") +include("prescribed_ocean.jl") + +# Grab the ocean state (Needs extending for new ocean models) +ocean_state(::Nothing) = (u = nothing, v = nothing, T = nothing, S = nothing) +ocean_state(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = ocean_state(ocean.model) +ocean_state(ocean) = (u = ocean.velocities.u, + v = ocean.velocities.v, + T = ocean.tracers.T, + S = ocean.tracers.S) end # module diff --git a/src/OceanSimulations/prescribed_ocean_model.jl b/src/OceanSimulations/prescribed_ocean.jl similarity index 83% rename from src/OceanSimulations/prescribed_ocean_model.jl rename to src/OceanSimulations/prescribed_ocean.jl index e1ff82296..661c77e6b 100644 --- a/src/OceanSimulations/prescribed_ocean_model.jl +++ b/src/OceanSimulations/prescribed_ocean.jl @@ -12,7 +12,7 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} +struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} architecture :: Arch grid :: G clock :: Clock{C} @@ -22,7 +22,7 @@ struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} end """ - PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` on a `grid` with a `clock`. @@ -34,7 +34,7 @@ Arguments following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. """ -function PrescribedOceanModel(timeseries=NamedTuple(); +function PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) @@ -60,14 +60,14 @@ function PrescribedOceanModel(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) end ##### ##### Need to extend a couple of methods ##### -function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) +function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) tick!(model.clock, Δt) time = Time(model.clock.time) @@ -99,8 +99,8 @@ function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) return nothing end -update_state!(::PrescribedOceanModel) = nothing -timestepper(::PrescribedOceanModel) = nothing +update_state!(::PrescribedOcean) = nothing +timestepper(::PrescribedOcean) = nothing -reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6 +reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 From 45d21912902f89ef379585108b2349778533f03c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:57:01 +0100 Subject: [PATCH 49/56] back how it was for the moment --- .../atmosphere_ocean_fluxes.jl | 5 ++++- src/OceanSimulations/OceanSimulations.jl | 12 ++---------- ...ibed_ocean.jl => prescribed_ocean_model.jl} | 18 +++++++++--------- 3 files changed, 15 insertions(+), 20 deletions(-) rename src/OceanSimulations/{prescribed_ocean.jl => prescribed_ocean_model.jl} (83%) diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl index 4a1b2c830..eac162d78 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl @@ -12,7 +12,10 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) arch = architecture(grid) clock = coupled_model.clock - ocean_state = ocean_state(ocean) + ocean_state = (u = ocean.model.u, + v = ocean.model.v, + T = ocean.model.T, + S = ocean.model.S) atmosphere_fields = coupled_model.interfaces.exchanger.exchange_atmosphere_state diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 353160a04..49d48aa6e 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,6 @@ module OceanSimulations -export ocean_simulation, PrescribedOcean, ocean_state +export ocean_simulation, PrescribedOceanModel using Oceananigans using Oceananigans.Units @@ -50,14 +50,6 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("ocean_simulation.jl") -include("prescribed_ocean.jl") - -# Grab the ocean state (Needs extending for new ocean models) -ocean_state(::Nothing) = (u = nothing, v = nothing, T = nothing, S = nothing) -ocean_state(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = ocean_state(ocean.model) -ocean_state(ocean) = (u = ocean.velocities.u, - v = ocean.velocities.v, - T = ocean.tracers.T, - S = ocean.tracers.S) +include("prescribed_ocean_model.jl") end # module diff --git a/src/OceanSimulations/prescribed_ocean.jl b/src/OceanSimulations/prescribed_ocean_model.jl similarity index 83% rename from src/OceanSimulations/prescribed_ocean.jl rename to src/OceanSimulations/prescribed_ocean_model.jl index 661c77e6b..e1ff82296 100644 --- a/src/OceanSimulations/prescribed_ocean.jl +++ b/src/OceanSimulations/prescribed_ocean_model.jl @@ -12,7 +12,7 @@ import Oceananigans.Architectures: on_architecture ##### A prescribed ocean... ##### -struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} +struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} architecture :: Arch grid :: G clock :: Clock{C} @@ -22,7 +22,7 @@ struct PrescribedOcean{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch} end """ - PrescribedOcean(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) + PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel` on a `grid` with a `clock`. @@ -34,7 +34,7 @@ Arguments following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries` and reside on the provided `grid`. """ -function PrescribedOcean(timeseries=NamedTuple(); +function PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0)) @@ -60,14 +60,14 @@ function PrescribedOcean(timeseries=NamedTuple(); T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ))) S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ))) - return PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) + return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries) end ##### ##### Need to extend a couple of methods ##### -function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) +function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true) tick!(model.clock, Δt) time = Time(model.clock.time) @@ -99,8 +99,8 @@ function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) return nothing end -update_state!(::PrescribedOcean) = nothing -timestepper(::PrescribedOcean) = nothing +update_state!(::PrescribedOceanModel) = nothing +timestepper(::PrescribedOceanModel) = nothing -reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6 -heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6 +reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6 +heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6 From 9eff489db63b0091df68c71085a71e39e67d3748 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 15:59:23 +0100 Subject: [PATCH 50/56] see how it goes --- .../InterfaceComputations/atmosphere_ocean_fluxes.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl index eac162d78..89b736db2 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl @@ -12,10 +12,10 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) arch = architecture(grid) clock = coupled_model.clock - ocean_state = (u = ocean.model.u, - v = ocean.model.v, - T = ocean.model.T, - S = ocean.model.S) + ocean_state = (u = ocean.model.velocities.u, + v = ocean.model.velocities.v, + T = ocean.model.tracers.T, + S = ocean.model.tracers.S) atmosphere_fields = coupled_model.interfaces.exchanger.exchange_atmosphere_state From f1c8bdead57ec85c43c6540bae44f26a7886c72f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 17:26:40 +0100 Subject: [PATCH 51/56] just add this for the moment --- .../InterfaceComputations/assemble_net_fluxes.jl | 16 +++++++++++----- .../atmosphere_ocean_fluxes.jl | 1 - .../atmosphere_sea_ice_fluxes.jl | 2 +- .../component_interfaces.jl | 2 +- 4 files changed, 13 insertions(+), 8 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 37bbb29c9..e415628d7 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -173,7 +173,8 @@ function compute_net_sea_ice_fluxes!(coupled_model) kernel_parameters = interface_kernel_parameters(grid) - sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_ocean_interface.temperature + sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature + ice_concentration = sea_ice_concentration(sea_ice) launch!(arch, grid, kernel_parameters, _assemble_net_sea_ice_fluxes!, @@ -186,6 +187,7 @@ function compute_net_sea_ice_fluxes!(coupled_model) freshwater_flux, sea_ice_surface_temperature, downwelling_radiation, + ice_concentration, sea_ice_properties, atmos_sea_ice_properties) @@ -201,6 +203,7 @@ end freshwater_flux, # Where do we add this one? surface_temperature, downwelling_radiation, + ice_concentration, sea_ice_properties, atmos_sea_ice_properties) @@ -211,13 +214,13 @@ end @inbounds begin Ts = surface_temperature[i, j, kᴺ] Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts) - + ℵi = ice_concentration[i, j, 1] Qs = downwelling_radiation.Qs[i, j, 1] Qℓ = downwelling_radiation.Qℓ[i, j, 1] Qc = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux Qv = atmosphere_sea_ice_fluxes.latent_heat[i, j, 1] # latent heat flux Qf = sea_ice_ocean_fluxes.frazil_heat[i, j, 1] # frazil heat flux - Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1] # interfacial heat flux + Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1] # interfacial heat flux end ρτx = atmosphere_sea_ice_fluxes.x_momentum # zonal momentum flux @@ -230,12 +233,15 @@ end Qu = upwelling_radiation(Ts, σ, ϵ) Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ) - ΣQt = Qd + Qu + Qc + Qv - ΣQb = Qf + Qi + ΣQt = (Qd + Qu + Qc + Qv) #* ℵi + ΣQb = Qi # Mask fluxes over land for convenience inactive = inactive_node(i, j, kᴺ, grid, c, c, c) + if ΣQt - Qc - Qv - Qu > 0 + @show ΣQt, Qc, Qv, Qu, Qd + end @inbounds top_fluxes.heat[i, j, 1] = ifelse(inactive, zero(grid), ΣQt) @inbounds top_fluxes.u[i, j, 1] = ifelse(inactive, zero(grid), ℑxᶠᵃᵃ(i, j, 1, grid, ρτx)) @inbounds top_fluxes.v[i, j, 1] = ifelse(inactive, zero(grid), ℑyᵃᶠᵃ(i, j, 1, grid, ρτy)) diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl index 89b736db2..aff140335 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl @@ -71,7 +71,6 @@ end i, j = @index(Global, NTuple) kᴺ = size(grid, 3) # index of the top ocean cell - time = Time(clock.time) @inbounds begin uₐ = atmosphere_state.u[i, j, 1] diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index c09e60f32..0bc5f1c72 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -116,7 +116,7 @@ end h_bℓ = atmosphere_state.h_bℓ) downwelling_radiation = (; Qs, Qℓ) - local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ) + local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ, ℵ=ℵᵢ) # Estimate initial interface state u★ = convert(FT, 1e-4) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 9b0f49fd4..cdacf9cca 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -266,7 +266,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)), - atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)), + atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid), stability_functions=atmosphere_sea_ice_stability_functions()), atmosphere_ocean_interface_temperature = BulkTemperature(), atmosphere_ocean_velocity_difference = RelativeVelocity(), atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean), From e1ae3040dfadc8b71e8c655365ccd4799ee0253f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 17:34:32 +0100 Subject: [PATCH 52/56] we need a looot of iterations for the moment --- .../InterfaceComputations/assemble_net_fluxes.jl | 7 ++----- .../InterfaceComputations/component_interfaces.jl | 2 +- .../InterfaceComputations/interface_states.jl | 2 +- 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index e415628d7..9e7a174bc 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -233,15 +233,12 @@ end Qu = upwelling_radiation(Ts, σ, ϵ) Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ) - ΣQt = (Qd + Qu + Qc + Qv) #* ℵi - ΣQb = Qi + ΣQt = (Qd + Qu + Qc + Qv) # * ℵi We need to multiply these times the concentration? + ΣQb = Qf + Qi # Mask fluxes over land for convenience inactive = inactive_node(i, j, kᴺ, grid, c, c, c) - if ΣQt - Qc - Qv - Qu > 0 - @show ΣQt, Qc, Qv, Qu, Qd - end @inbounds top_fluxes.heat[i, j, 1] = ifelse(inactive, zero(grid), ΣQt) @inbounds top_fluxes.u[i, j, 1] = ifelse(inactive, zero(grid), ℑxᶠᵃᵃ(i, j, 1, grid, ρτx)) @inbounds top_fluxes.v[i, j, 1] = ifelse(inactive, zero(grid), ℑyᵃᶠᵃ(i, j, 1, grid, ρτy)) diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index cdacf9cca..0336dd3d0 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -266,7 +266,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice; radiation = Radiation(), freshwater_density = 1000, atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)), - atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid), stability_functions=atmosphere_sea_ice_stability_functions()), + atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid), solver_maxiter=10000, stability_functions=atmosphere_sea_ice_stability_functions()), atmosphere_ocean_interface_temperature = BulkTemperature(), atmosphere_ocean_velocity_difference = RelativeVelocity(), atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean), diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 108c9413e..5aa4d25c8 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -225,7 +225,7 @@ end Tₛ = (Tᵢ - h / k * (Qₐ + 4α * Tₛ⁻^4)) / (1 + 4α * h * Tₛ⁻^3 / k) =# - T★ = Tᵢ - Qₐ * h / k + T★ = Tᵢ - Qₐ * h / k # * ℵi We need to multiply the flux by the concentration? # Fix a NaN T★ = ifelse(isnan(T★), Tₛ⁻, T★) From af3a2e73268655b600247bfd1dcc870384e6e4ae Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 13 Mar 2025 17:54:43 +0100 Subject: [PATCH 53/56] try it like this --- .../InterfaceComputations/assemble_net_fluxes.jl | 2 +- src/OceanSeaIceModels/InterfaceComputations/interface_states.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 9e7a174bc..7391dfd70 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -233,7 +233,7 @@ end Qu = upwelling_radiation(Ts, σ, ϵ) Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ) - ΣQt = (Qd + Qu + Qc + Qv) # * ℵi We need to multiply these times the concentration? + ΣQt = (Qd + Qu + Qc + Qv) * ℵi # We need to multiply these times the concentration? ΣQb = Qf + Qi # Mask fluxes over land for convenience diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 5aa4d25c8..d83713e4a 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -225,7 +225,7 @@ end Tₛ = (Tᵢ - h / k * (Qₐ + 4α * Tₛ⁻^4)) / (1 + 4α * h * Tₛ⁻^3 / k) =# - T★ = Tᵢ - Qₐ * h / k # * ℵi We need to multiply the flux by the concentration? + T★ = Tᵢ - Qₐ * h / k * ℵi # We need to multiply the flux by the concentration? # Fix a NaN T★ = ifelse(isnan(T★), Tₛ⁻, T★) From e89e9f6babf06057ff4dae0c4b9011874baab7b6 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Thu, 13 Mar 2025 13:10:53 -0400 Subject: [PATCH 54/56] small change --- experiments/arctic_simulation.jl | 2 +- src/OceanSeaIceModels/ocean_sea_ice_model.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 6538434c9..449affd31 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -81,7 +81,7 @@ dynamics = SeaIceMomentumEquation(grid; rheology = ElastoViscoPlasticRheology(), solver = SplitExplicitSolver(120)) -sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition) #, dynamics, advection=WENO(order=7)) set!(sea_ice.model, h=Metadata(:sea_ice_thickness; dataset), ℵ=Metadata(:sea_ice_concentration; dataset)) diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 418ec41fe..81386b13f 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -164,7 +164,7 @@ function above_freezing_ocean_temperature!(ocean, sea_ice::SeaIceSimulation) grid = ocean.model.grid arch = architecture(grid) - launch!(arch, grid, :xy, _above_freezing_ocean_temperature!, T, grid, S, ℵ, liquidus) + launch!(arch, grid, :xy, _above_freezing_ocean_temperature!, T, grid, S, liquidus) return nothing end From bc74c087e2ca89727eb7934ae9ee4c2f2658e662 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Thu, 13 Mar 2025 13:19:55 -0400 Subject: [PATCH 55/56] bugfix --- experiments/arctic_simulation.jl | 2 +- .../InterfaceComputations/interface_states.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 449affd31..13bbb7224 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -12,7 +12,7 @@ using Printf using CUDA CUDA.device!(1) -arch = GPU() +arch = CPU() r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=2000) z_faces = MutableVerticalDiscretization(r_faces) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index d83713e4a..8a693b54f 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -212,12 +212,13 @@ end F = st.internal_flux k = F.conductivity h = Ψᵢ.h + ℵ = Ψᵢ.ℵ # Bottom temperature at the melting temperature Tᵢ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) Tᵢ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ) Tₛ⁻ = Ψₛ.T - + #= σ = ℙₛ.radiation.σ ϵ = ℙₛ.radiation.ϵ @@ -225,7 +226,7 @@ end Tₛ = (Tᵢ - h / k * (Qₐ + 4α * Tₛ⁻^4)) / (1 + 4α * h * Tₛ⁻^3 / k) =# - T★ = Tᵢ - Qₐ * h / k * ℵi # We need to multiply the flux by the concentration? + T★ = Tᵢ - Qₐ * h / k * ℵ # We need to multiply the flux by the concentration? # Fix a NaN T★ = ifelse(isnan(T★), Tₛ⁻, T★) From 57508b39a3cd108732f8b47e08b1b40dc87b7244 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Thu, 13 Mar 2025 13:27:44 -0400 Subject: [PATCH 56/56] back to gpu --- experiments/arctic_simulation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 13bbb7224..449affd31 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -12,7 +12,7 @@ using Printf using CUDA CUDA.device!(1) -arch = CPU() +arch = GPU() r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=2000) z_faces = MutableVerticalDiscretization(r_faces)