From 1667e33916eef85aea7505a8c7e7df86ea8ca523 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 30 Apr 2025 16:28:51 +0200 Subject: [PATCH 01/19] start climatoloy --- .../flux_climatology/flux_climatology.jl | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 897fce64..c0a97680 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -1,6 +1,5 @@ using ClimaOcean using ClimaOcean.ECCO -using ClimaOcean.ECCO: all_ECCO_dates using Oceananigans using Oceananigans.Utils using Oceananigans.Fields: ZeroField, location, VelocityFields @@ -142,15 +141,21 @@ function PrescribedOcean(timeseries; end # ...with prescribed velocity and tracer fields -version = ECCO4Monthly() -dates = all_ECCO_dates(version)[1:24] +dataset = ECCO4Monthly() +arch = CPU() -arch = CPU() +start_date = DateTime(1992, 1, 1) +end_date = DateTime(1993, 1, 1) -u = ECCOFieldTimeSeries(:u_velocity, version; time_indices_in_memory = 13, architecture = arch, dates) -v = ECCOFieldTimeSeries(:v_velocity, version; time_indices_in_memory = 13, architecture = arch, dates) -T = ECCOFieldTimeSeries(:temperature, version; time_indices_in_memory = 13, architecture = arch, dates) -S = ECCOFieldTimeSeries(:salinity, version; time_indices_in_memory = 13, architecture = arch, dates) +T_meta = Metadata(:temperature; start_date, end_date, dataset) +S_meta = Metadata(:salinity; start_date, end_date, dataset) +u_meta = Metadata(:u_velocity; start_date, end_date, dataset) +v_meta = Metadata(:v_velocity; start_date, end_date, dataset) + +u = FieldTimeSeries(u_meta, arch; time_indices_in_memory = 13) +v = FieldTimeSeries(v_meta, arch; time_indices_in_memory = 13) +T = FieldTimeSeries(T_meta, arch; time_indices_in_memory = 13) +S = FieldTimeSeries(S_meta, arch; time_indices_in_memory = 13) grid = ECCO.ECCO_immersed_grid(arch) From e26b67e906caceaa3fe0def3a296a3f0e2baf64f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 30 Apr 2025 16:44:51 +0200 Subject: [PATCH 02/19] few fixes + sensible and latent heat fluxes --- experiments/flux_climatology/flux_climatology.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index c0a97680..d2da05bf 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -1,5 +1,6 @@ using ClimaOcean using ClimaOcean.ECCO +using Dates using Oceananigans using Oceananigans.Utils using Oceananigans.Fields: ZeroField, location, VelocityFields @@ -84,7 +85,11 @@ function compute_flux_climatology(earth) Jᵀ = FluxStatistics(net_fluxes.T) Jˢ = FluxStatistics(net_fluxes.S) - stats = (; τx, τy, Jᵀ, Jˢ) + atmos_ocean_fluxes = earth.model.interfaces.atmosphere_ocean_interface.fluxes + Qc = FluxStatistics(atmos_ocean_fluxes.sensible_heat) + Qv = FluxStatistics(atmos_ocean_fluxes.latent_heat) + + stats = (; τx, τy, Jᵀ, Jˢ, Qc, Qv) function update_flux_stats!(earth) Δt = earth.Δt @@ -114,7 +119,7 @@ end ##### A prescribed ocean... ##### -struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing} +struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing, A} architecture :: A grid :: G clock :: Clock{C} @@ -232,5 +237,3 @@ end add_callback!(earth, progress, IterationInterval(10)) stats = compute_flux_climatology(earth) - - From 63bd72a79b841a3b7166f1915b66dfa04b34ec32 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 1 May 2025 10:11:22 +0200 Subject: [PATCH 03/19] no space --- experiments/flux_climatology/flux_climatology.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index d2da05bf..5a1ff9c0 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -157,10 +157,10 @@ S_meta = Metadata(:salinity; start_date, end_date, dataset) u_meta = Metadata(:u_velocity; start_date, end_date, dataset) v_meta = Metadata(:v_velocity; start_date, end_date, dataset) -u = FieldTimeSeries(u_meta, arch; time_indices_in_memory = 13) -v = FieldTimeSeries(v_meta, arch; time_indices_in_memory = 13) -T = FieldTimeSeries(T_meta, arch; time_indices_in_memory = 13) -S = FieldTimeSeries(S_meta, arch; time_indices_in_memory = 13) +u = FieldTimeSeries(u_meta, arch; time_indices_in_memory=13) +v = FieldTimeSeries(v_meta, arch; time_indices_in_memory=13) +T = FieldTimeSeries(T_meta, arch; time_indices_in_memory=13) +S = FieldTimeSeries(S_meta, arch; time_indices_in_memory=13) grid = ECCO.ECCO_immersed_grid(arch) From 9760ef85306536ffb6b5d9d31f687d4780c7e864 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 1 May 2025 11:21:12 +0200 Subject: [PATCH 04/19] some fixing --- experiments/flux_climatology/flux_climatology.jl | 2 +- .../InterfaceComputations/interpolate_atmospheric_state.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 5a1ff9c0..1b1b7af3 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -165,7 +165,7 @@ S = FieldTimeSeries(S_meta, arch; time_indices_in_memory=13) grid = ECCO.ECCO_immersed_grid(arch) ocean_model = PrescribedOcean((; u, v, T, S); grid) -ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days) +ocean = Simulation(ocean_model, Δt=3hours, stop_time=365days) ##### ##### Need to extend a couple of methods diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index 6997df68..28d0d5f4 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -111,8 +111,8 @@ function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosph # # TODO: find a better design for this that doesn't have redundant # arrays for the barotropic potential - u_potential = forcing_barotropic_potential(ocean.model.forcing.u) - v_potential = forcing_barotropic_potential(ocean.model.forcing.v) + u_potential = forcing_barotropic_potential(ocean) + v_potential = forcing_barotropic_potential(ocean) ρₒ = coupled_model.interfaces.ocean_properties.reference_density if !isnothing(u_potential) From c8f27ad1519841fd3c64fdbfe1f981edd062302a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 1 May 2025 11:45:39 +0200 Subject: [PATCH 05/19] run it on GPU for the moment --- experiments/flux_climatology/flux_climatology.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 1b1b7af3..2d8691ca 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -147,7 +147,7 @@ end # ...with prescribed velocity and tracer fields dataset = ECCO4Monthly() -arch = CPU() +arch = GPU() start_date = DateTime(1992, 1, 1) end_date = DateTime(1993, 1, 1) @@ -173,7 +173,7 @@ ocean = Simulation(ocean_model, Δt=3hours, stop_time=365days) function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true) tick!(model.clock, Δt) - time = Time(model.clock.time) + time = Oceananigans.Units.Time(model.clock.time) possible_fts = merge(model.velocities, model.tracers) time_series_tuple = extract_field_time_series(possible_fts) From 5edabc084acfce7dc02050382ceeb86d847e89d9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 May 2025 14:16:13 +0300 Subject: [PATCH 06/19] add dataset_type --- src/DataWrangling/DataWrangling.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/DataWrangling/DataWrangling.jl b/src/DataWrangling/DataWrangling.jl index 1b27806e..7c413d74 100644 --- a/src/DataWrangling/DataWrangling.jl +++ b/src/DataWrangling/DataWrangling.jl @@ -20,6 +20,8 @@ using KernelAbstractions: @kernel, @index using Oceananigans.DistributedComputations using Adapt +const dataset_defaults = Oceananigans.Defaults(FloatType=Float32) + ##### ##### Downloading utilities ##### From 9ccb0cc90d4ab16e3ff54b1e0bd9eac3c006d5b5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 May 2025 14:16:32 +0300 Subject: [PATCH 07/19] use dataset_defaults.FloatType --- src/DataWrangling/metadata.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/DataWrangling/metadata.jl b/src/DataWrangling/metadata.jl index 13c25470..6b400975 100644 --- a/src/DataWrangling/metadata.jl +++ b/src/DataWrangling/metadata.jl @@ -227,7 +227,7 @@ function dataset_latitude_extent end Return an empty field of `metadata` on `architecture` and with `horizontal_halo`s. """ -function empty_field(metadata::Metadata; +function empty_field(metadata::Metadata, FT=dataset_defaults.FloatType; architecture = CPU(), horizontal_halo = (7, 7)) @@ -251,7 +251,7 @@ function empty_field(metadata::Metadata; sz = (Nx, Ny) end - grid = LatitudeLongitudeGrid(architecture, Float32; halo, longitude, latitude, z, + grid = LatitudeLongitudeGrid(architecture, FT; halo, longitude, latitude, z, size = sz, topology = (TX, TY, TZ)) return Field{loc...}(grid) From 894ab394cd88e1a26b04fdd08678595ced1c2bbb Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 May 2025 14:16:54 +0300 Subject: [PATCH 08/19] add Qc+Qv + better averaging that do not depend on the end_time --- .../flux_climatology/flux_climatology.jl | 135 +++++++++++------- .../flux_climatology/visualize_climatology.jl | 111 ++++++++------ 2 files changed, 144 insertions(+), 102 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 2d8691ca..c1fce730 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -19,61 +19,72 @@ import Oceananigans.Models: timestepper, update_model_field_time_series! import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity import Oceananigans.Architectures: on_architecture +ClimaOcean.DataWrangling.dataset_defaults.FloatType = Float64 + ##### ##### A Data structure that holds flux statistics ##### struct FluxStatistics{F} - avg :: F + mean :: F + meansq :: F std :: F max :: F min :: F end function FluxStatistics(f::Field) - a = similar(f) - s = similar(f) - p = similar(f) - m = similar(f) - - fill!(a, 0) - fill!(s, 0) - fill!(p, 0) - fill!(m, 0) - - return FluxStatistics(a, s, p, m) + mean = similar(f) + meansq = similar(f) + std = similar(f) + max = similar(f) + min = similar(f) + + fill!(mean, 0) + fill!(meansq, 0) + fill!(std, 0) + fill!(max, 0) + fill!(min, 0) + + return FluxStatistics(mean, meansq, std, max, min) end -Adapt.adapt_structure(to, f::FluxStatistics) = FluxStatistics(Adapt.adapt(to, f.avg), +compute_std!(f::FluxStatistics) = parent(f.std) .= sqrt.(parent(f.meansq) .- parent(f.mean).^2) + +Adapt.adapt_structure(to, f::FluxStatistics) = FluxStatistics(Adapt.adapt(to, f.mean), + Adapt.adapt(to, f.meansq), Adapt.adapt(to, f.std), Adapt.adapt(to, f.max), Adapt.adapt(to, f.min)) -on_architecture(arch, f::FluxStatistics) = FluxStatistics(on_architecture(arch, f.avg), +on_architecture(arch, f::FluxStatistics) = FluxStatistics(on_architecture(arch, f.mean), + on_architecture(arch, f.meansq), on_architecture(arch, f.std), on_architecture(arch, f.max), on_architecture(arch, f.min)) -function update_stats!(stats::FluxStatistics, flux, Δt, stop_time) +function update_stats!(stats::FluxStatistics, flux, iteration) grid = flux.grid arch = architecture(grid) - launch!(arch, grid, :xy, _update_stats!, stats, flux, Δt, stop_time) + launch!(arch, grid, :xy, _update_stats!, stats, flux, iteration) end -@kernel function _update_stats!(stats, flux, Δt, stop_time) +@kernel function _update_stats!(stats, flux, iteration) i, j = @index(Global, NTuple) + inverse_iteration = 1 / (iteration + 1) + @inbounds begin - f = flux[i, j, 1] - stats.avg[i, j, 1] += f * Δt / stop_time - stats.std[i, j, 1] += f^2 * Δt / stop_time + f = flux[i, j, 1] + stats.mean[i, j, 1] *= 1 - inverse_iteration + stats.mean[i, j, 1] += f * inverse_iteration + stats.meansq[i, j, 1] *= 1 - inverse_iteration + stats.meansq[i, j, 1] += f^2 * inverse_iteration stats.max[i, j, 1] = max(stats.max[i, j, 1], f) stats.min[i, j, 1] = min(stats.min[i, j, 1], f) end end -finalize_std!(f::FluxStatistics) = parent(f.std) .= sqrt.(parent(f.std) .- parent(f.avg).^2) - ##### ##### A function to compute flux statistics ##### @@ -92,13 +103,13 @@ function compute_flux_climatology(earth) stats = (; τx, τy, Jᵀ, Jˢ, Qc, Qv) function update_flux_stats!(earth) - Δt = earth.Δt - stop_time = earth.stop_time - - update_stats!(τx, net_fluxes.u, Δt, stop_time) - update_stats!(τy, net_fluxes.v, Δt, stop_time) - update_stats!(Jᵀ, net_fluxes.T, Δt, stop_time) - update_stats!(Jˢ, net_fluxes.S, Δt, stop_time) + iteration = earth.model.clock.iteration + update_stats!(τx, net_fluxes.u, iteration) + update_stats!(τy, net_fluxes.v, iteration) + update_stats!(Jᵀ, net_fluxes.T, iteration) + update_stats!(Jˢ, net_fluxes.S, iteration) + update_stats!(Qc, atmos_ocean_fluxes.sensible_heat, iteration) + update_stats!(Qv, atmos_ocean_fluxes.latent_heat, iteration) return nothing end @@ -107,10 +118,12 @@ function compute_flux_climatology(earth) run!(earth) - finalize_std!(τx) - finalize_std!(τy) - finalize_std!(Jᵀ) - finalize_std!(Jˢ) + compute_std!(τx) + compute_std!(τy) + compute_std!(Jᵀ) + compute_std!(Jˢ) + update_stats!(Qc) + update_stats!(Qv) return stats end @@ -147,25 +160,29 @@ end # ...with prescribed velocity and tracer fields dataset = ECCO4Monthly() -arch = GPU() +arch = CPU() start_date = DateTime(1992, 1, 1) -end_date = DateTime(1993, 1, 1) +end_date = DateTime(1992, 12, 31) + +stop_time = Day(end_date - start_date).value * Oceananigans.Units.days + +time_indices_in_memory = 13 -T_meta = Metadata(:temperature; start_date, end_date, dataset) -S_meta = Metadata(:salinity; start_date, end_date, dataset) u_meta = Metadata(:u_velocity; start_date, end_date, dataset) v_meta = Metadata(:v_velocity; start_date, end_date, dataset) +T_meta = Metadata(:temperature; start_date, end_date, dataset) +S_meta = Metadata(:salinity; start_date, end_date, dataset) -u = FieldTimeSeries(u_meta, arch; time_indices_in_memory=13) -v = FieldTimeSeries(v_meta, arch; time_indices_in_memory=13) -T = FieldTimeSeries(T_meta, arch; time_indices_in_memory=13) -S = FieldTimeSeries(S_meta, arch; time_indices_in_memory=13) +u = FieldTimeSeries(u_meta, arch; time_indices_in_memory) +v = FieldTimeSeries(v_meta, arch; time_indices_in_memory) +T = FieldTimeSeries(T_meta, arch; time_indices_in_memory) +S = FieldTimeSeries(S_meta, arch; time_indices_in_memory) grid = ECCO.ECCO_immersed_grid(arch) ocean_model = PrescribedOcean((; u, v, T, S); grid) -ocean = Simulation(ocean_model, Δt=3hours, stop_time=365days) +ocean = Simulation(ocean_model; Δt=3hours, stop_time) ##### ##### Need to extend a couple of methods @@ -206,34 +223,42 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(1000)) ##### A prescribed earth... ##### -earth_model = OceanSeaIceModel(ocean, nothing; atmosphere, radiation = Radiation(ocean_emissivity=0, ocean_albedo=1)) -earth = Simulation(earth_model, Δt=3hours, stop_time=365days) +earth_model = OceanSeaIceModel(ocean, nothing; atmosphere, radiation = Radiation(arch)) +earth = Simulation(earth_model; Δt=3hours, stop_time) wall_time = Ref(time_ns()) # Just checking that the ocean evolves as expected function progress(sim) - ocean = sim.model.ocean - u, v, w = ocean.model.velocities - T = ocean.model.tracers.T + net_fluxes = sim.model.interfaces.net_fluxes.ocean_surface + atmos_ocean_fluxes = sim.model.interfaces.atmosphere_ocean_interface.fluxes + + Q = net_fluxes.Q + Qc = atmos_ocean_fluxes.sensible_heat + Qv = atmos_ocean_fluxes.latent_heat - Tmax = maximum(interior(T)) - Tmin = minimum(interior(T)) + τx = net_fluxes.u + τy = net_fluxes.v - umax = (maximum(abs, interior(u)), - maximum(abs, interior(v))) + Qmax = maximum(Q) + Qmin = minimum(Q) + Qcmax = maximum(Qc) + Qcmin = minimum(Qc) + Qvmax = maximum(Qv) + Qvmin = minimum(Qv) + τmax = (maximum(abs, τx), maximum(abs, τy)) step_time = 1e-9 * (time_ns() - wall_time[]) - msg = @sprintf("Iter: %d, time: %s, Δt: %s", iteration(sim), prettytime(sim), prettytime(sim.Δt)) - msg *= @sprintf(", max|u|: (%.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s", - umax..., Tmax, Tmin, prettytime(step_time)) + msg = @sprintf("Iter: %d, time: %s", iteration(sim), prettytime(sim)) + msg *= @sprintf(", max|τ|: (%.2e, %.2e) m² s⁻², extrema(Q): (%.2f, %.2f) W m⁻², extrema(Qc): (%.2f, %.2f) W m⁻², extrema(Qv): (%.2f, %.2f) W m⁻², wall time: %s", + τmax..., Qmin, Qmax, Qcmin, Qcmax, Qvmin, Qvmax, prettytime(step_time)) @info msg wall_time[] = time_ns() end -add_callback!(earth, progress, IterationInterval(10)) +add_callback!(earth, progress, IterationInterval(16)) stats = compute_flux_climatology(earth) diff --git a/experiments/flux_climatology/visualize_climatology.jl b/experiments/flux_climatology/visualize_climatology.jl index 7da86c0c..4382e057 100644 --- a/experiments/flux_climatology/visualize_climatology.jl +++ b/experiments/flux_climatology/visualize_climatology.jl @@ -1,86 +1,103 @@ - -using CairoMakie - -fig = Figure(size = (1200, 600)) -axQ = Axis(fig[1, 1], title = "Heat Flux") -axS = Axis(fig[1, 2], title = "Salinity Flux") -axu = Axis(fig[2, 1], title = "Zonal stress") -axv = Axis(fig[2, 2], title = "Meridional stress") - ρ = reference_density(earth.model.ocean) cp = heat_capacity(earth.model.ocean) -hmQ = heatmap!(axQ, stats.Jᵀ.avg * ρ * cp, colormap = :balance, colorrange = (-100, 250)) -hmS = heatmap!(axS, stats.Jˢ.avg, colormap = :balance, colorrange = (-2.5e-6, 2.5e-6)) -hmu = heatmap!(axu, stats.τx.avg * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) -hmv = heatmap!(axv, stats.τy.avg * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) +using GLMakie -Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") -Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") -Colorbar(fig[2, 0], hmu, flipaxis=false, label = "N/m²") -Colorbar(fig[2, 3], hmv, flipaxis=true, label = "N/m²") +fig = Figure(size = (1200, 900)) +axQ = Axis(fig[1, 1], title = "Heat Flux") +axS = Axis(fig[1, 2], title = "Salinity Flux") +axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") +axQv = Axis(fig[2, 2], title = "Latent Heat Flux") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") + +hmQ = heatmap!(axQ, stats.Jᵀ.mean * ρ * cp, colormap = :balance, colorrange = (-100, 250)) +hmS = heatmap!(axS, stats.Jˢ.mean, colormap = :balance, colorrange = (-2.5e-6, 2.5e-6)) +hmQc = heatmap!(axQc, stats.Qc.mean, colormap = :balance, colorrange = (-100, 250)) +hmQv = heatmap!(axQv, stats.Qv.mean, colormap = :balance, colorrange = (-100, 250)) +hmu = heatmap!(axu, stats.τx.mean * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) +hmv = heatmap!(axv, stats.τy.mean * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) + +Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") +Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") +Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") +Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") save("mean_fluxes.png", fig) -fig = Figure(size = (1200, 600)) + +fig = Figure(size = (1200, 900)) axQ = Axis(fig[1, 1], title = "Heat Flux") axS = Axis(fig[1, 2], title = "Salinity Flux") -axu = Axis(fig[2, 1], title = "Zonal stress") -axv = Axis(fig[2, 2], title = "Meridional stress") - -ρ = reference_density(earth.model.ocean) -cp = heat_capacity(earth.model.ocean) +axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") +axQv = Axis(fig[2, 2], title = "Latent Heat Flux") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") hmQ = heatmap!(axQ, stats.Jᵀ.std * ρ * cp, colormap = :deep, colorrange = (0, 200)) hmS = heatmap!(axS, stats.Jˢ.std, colormap = :deep, colorrange = (0, 5e-6)) +hmQc = heatmap!(axQc, stats.Qc.std, colormap = :deep, colorrange = (0, 200)) +hmQv = heatmap!(axQv, stats.Qv.std, colormap = :deep, colorrange = (0, 200)) hmu = heatmap!(axu, stats.τx.std * ρ, colormap = :deep, colorrange = (0, 0.3)) hmv = heatmap!(axv, stats.τy.std * ρ, colormap = :deep, colorrange = (0, 0.3)) -Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") -Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") -Colorbar(fig[2, 0], hmu, flipaxis=false, label = "N/m²") -Colorbar(fig[2, 3], hmv, flipaxis=true, label = "N/m²") +Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") +Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") +Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") +Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") save("fluxes_std.png", fig) -fig = Figure(size = (1200, 600)) + +fig = Figure(size = (1200, 900)) axQ = Axis(fig[1, 1], title = "Heat Flux") axS = Axis(fig[1, 2], title = "Salinity Flux") -axu = Axis(fig[2, 1], title = "Zonal stress") -axv = Axis(fig[2, 2], title = "Meridional stress") - -ρ = reference_density(earth.model.ocean) -cp = heat_capacity(earth.model.ocean) +axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") +axQv = Axis(fig[2, 2], title = "Latent Heat Flux") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") hmQ = heatmap!(axQ, stats.Jᵀ.max * ρ * cp, colormap = :magma, colorrange = (0, 1000)) hmS = heatmap!(axS, stats.Jˢ.max, colormap = :magma, colorrange = (1e-6, 5e-5)) +hmQc = heatmap!(axQc, stats.Qc.max, colormap = :magma, colorrange = (0, 1000)) +hmQv = heatmap!(axQv, stats.Qv.max, colormap = :magma, colorrange = (0, 1000)) hmu = heatmap!(axu, stats.τx.max * ρ, colormap = :magma, colorrange = (0.1, 2.0)) hmv = heatmap!(axv, stats.τy.max * ρ, colormap = :magma, colorrange = (0.1, 2.0)) -Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") -Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") -Colorbar(fig[2, 0], hmu, flipaxis=false, label = "N/m²") -Colorbar(fig[2, 3], hmv, flipaxis=true, label = "N/m²") +Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") +Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") +Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") +Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") save("fluxes_max.png", fig) -fig = Figure(size = (1200, 600)) + +fig = Figure(size = (1200, 900)) axQ = Axis(fig[1, 1], title = "Heat Flux") axS = Axis(fig[1, 2], title = "Salinity Flux") -axu = Axis(fig[2, 1], title = "Zonal stress") -axv = Axis(fig[2, 2], title = "Meridional stress") - -ρ = reference_density(earth.model.ocean) -cp = heat_capacity(earth.model.ocean) +axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") +axQv = Axis(fig[2, 2], title = "Latent Heat Flux") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") hmQ = heatmap!(axQ, stats.Jᵀ.min * ρ * cp, colormap = :haline, colorrange = (-1500, 0)) hmS = heatmap!(axS, stats.Jˢ.min, colormap = :haline, colorrange = (-5e-6, 0)) +hmQc = heatmap!(axQc, stats.Qc.min, colormap = :haline, colorrange = (-1000, 0)) +hmQv = heatmap!(axQv, stats.Qv.min, colormap = :haline, colorrange = (-1000, 0)) hmu = heatmap!(axu, stats.τx.min * ρ, colormap = :haline, colorrange = (-2.0, 0.1)) hmv = heatmap!(axv, stats.τy.min * ρ, colormap = :haline, colorrange = (-2.0, 0.1)) -Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") -Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") -Colorbar(fig[2, 0], hmu, flipaxis=false, label = "N/m²") -Colorbar(fig[2, 3], hmv, flipaxis=true, label = "N/m²") +Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") +Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") +Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") +Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") save("fluxes_min.png", fig) From f1398040d89901c2d1f0a0fc698f3c49b9695024 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 May 2025 14:53:54 +0300 Subject: [PATCH 09/19] fixes --- experiments/flux_climatology/flux_climatology.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index c1fce730..c70fc081 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -49,8 +49,6 @@ function FluxStatistics(f::Field) return FluxStatistics(mean, meansq, std, max, min) end -compute_std!(f::FluxStatistics) = parent(f.std) .= sqrt.(parent(f.meansq) .- parent(f.mean).^2) - Adapt.adapt_structure(to, f::FluxStatistics) = FluxStatistics(Adapt.adapt(to, f.mean), Adapt.adapt(to, f.meansq), Adapt.adapt(to, f.std), @@ -67,6 +65,10 @@ function update_stats!(stats::FluxStatistics, flux, iteration) grid = flux.grid arch = architecture(grid) launch!(arch, grid, :xy, _update_stats!, stats, flux, iteration) + + parent(stats.std) .= sqrt.(parent(stats.meansq) .- parent(stats.mean).^2) + + return nothing end @kernel function _update_stats!(stats, flux, iteration) @@ -74,6 +76,8 @@ end inverse_iteration = 1 / (iteration + 1) + # use iterative averaging via + # mean = mean * (1 - 1/iteration) + xn / iteration @inbounds begin f = flux[i, j, 1] stats.mean[i, j, 1] *= 1 - inverse_iteration @@ -118,13 +122,6 @@ function compute_flux_climatology(earth) run!(earth) - compute_std!(τx) - compute_std!(τy) - compute_std!(Jᵀ) - compute_std!(Jˢ) - update_stats!(Qc) - update_stats!(Qv) - return stats end From e6428dab9b0f261e9c058b9ac2be09c73e5ead33 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 May 2025 17:58:55 +0300 Subject: [PATCH 10/19] update figs --- .../flux_climatology/visualize_climatology.jl | 100 ++++++++++-------- 1 file changed, 53 insertions(+), 47 deletions(-) diff --git a/experiments/flux_climatology/visualize_climatology.jl b/experiments/flux_climatology/visualize_climatology.jl index 4382e057..788891b3 100644 --- a/experiments/flux_climatology/visualize_climatology.jl +++ b/experiments/flux_climatology/visualize_climatology.jl @@ -1,97 +1,103 @@ +px_per_unit = 4 + ρ = reference_density(earth.model.ocean) cp = heat_capacity(earth.model.ocean) using GLMakie fig = Figure(size = (1200, 900)) -axQ = Axis(fig[1, 1], title = "Heat Flux") -axS = Axis(fig[1, 2], title = "Salinity Flux") + +axQ = Axis(fig[1, 1], title = "Heat Flux") +axS = Axis(fig[1, 2], title = "Salinity Flux") axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") axQv = Axis(fig[2, 2], title = "Latent Heat Flux") -axu = Axis(fig[3, 1], title = "Zonal stress") -axv = Axis(fig[3, 2], title = "Meridional stress") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") -hmQ = heatmap!(axQ, stats.Jᵀ.mean * ρ * cp, colormap = :balance, colorrange = (-100, 250)) -hmS = heatmap!(axS, stats.Jˢ.mean, colormap = :balance, colorrange = (-2.5e-6, 2.5e-6)) -hmQc = heatmap!(axQc, stats.Qc.mean, colormap = :balance, colorrange = (-100, 250)) -hmQv = heatmap!(axQv, stats.Qv.mean, colormap = :balance, colorrange = (-100, 250)) -hmu = heatmap!(axu, stats.τx.mean * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) -hmv = heatmap!(axv, stats.τy.mean * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) +hmQ = heatmap!(axQ, stats.Jᵀ.mean * ρ * cp, colormap = :balance, colorrange = (-100, 250)) +hmS = heatmap!(axS, stats.Jˢ.mean, colormap = :balance, colorrange = (-2.5e-6, 2.5e-6)) +hmQc = heatmap!(axQc, stats.Qc.mean, colormap = :balance, colorrange = (-100, 250)) +hmQv = heatmap!(axQv, stats.Qv.mean, colormap = :balance, colorrange = (-100, 250)) +hmu = heatmap!(axu, stats.τx.mean * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) +hmv = heatmap!(axv, stats.τy.mean * ρ, colormap = :balance, colorrange = (-0.2, 0.2)) Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") -Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=true, label = "W/m²") Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") -save("mean_fluxes.png", fig) +save("fluxes_mean.png", fig; px_per_unit) fig = Figure(size = (1200, 900)) -axQ = Axis(fig[1, 1], title = "Heat Flux") -axS = Axis(fig[1, 2], title = "Salinity Flux") + +axQ = Axis(fig[1, 1], title = "Heat Flux") +axS = Axis(fig[1, 2], title = "Salinity Flux") axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") axQv = Axis(fig[2, 2], title = "Latent Heat Flux") -axu = Axis(fig[3, 1], title = "Zonal stress") -axv = Axis(fig[3, 2], title = "Meridional stress") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") -hmQ = heatmap!(axQ, stats.Jᵀ.std * ρ * cp, colormap = :deep, colorrange = (0, 200)) -hmS = heatmap!(axS, stats.Jˢ.std, colormap = :deep, colorrange = (0, 5e-6)) -hmQc = heatmap!(axQc, stats.Qc.std, colormap = :deep, colorrange = (0, 200)) -hmQv = heatmap!(axQv, stats.Qv.std, colormap = :deep, colorrange = (0, 200)) -hmu = heatmap!(axu, stats.τx.std * ρ, colormap = :deep, colorrange = (0, 0.3)) -hmv = heatmap!(axv, stats.τy.std * ρ, colormap = :deep, colorrange = (0, 0.3)) +hmQ = heatmap!(axQ, stats.Jᵀ.std * ρ * cp, colormap = :deep, colorrange = (0, 200)) +hmS = heatmap!(axS, stats.Jˢ.std, colormap = :deep, colorrange = (0, 5e-6)) +hmQc = heatmap!(axQc, stats.Qc.std, colormap = :deep, colorrange = (0, 200)) +hmQv = heatmap!(axQv, stats.Qv.std, colormap = :deep, colorrange = (0, 200)) +hmu = heatmap!(axu, stats.τx.std * ρ, colormap = :deep, colorrange = (0, 0.3)) +hmv = heatmap!(axv, stats.τy.std * ρ, colormap = :deep, colorrange = (0, 0.3)) Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") -Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=true, label = "W/m²") Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") -save("fluxes_std.png", fig) +save("fluxes_std.png", fig; px_per_unit) fig = Figure(size = (1200, 900)) -axQ = Axis(fig[1, 1], title = "Heat Flux") -axS = Axis(fig[1, 2], title = "Salinity Flux") + +axQ = Axis(fig[1, 1], title = "Heat Flux") +axS = Axis(fig[1, 2], title = "Salinity Flux") axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") axQv = Axis(fig[2, 2], title = "Latent Heat Flux") -axu = Axis(fig[3, 1], title = "Zonal stress") -axv = Axis(fig[3, 2], title = "Meridional stress") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") -hmQ = heatmap!(axQ, stats.Jᵀ.max * ρ * cp, colormap = :magma, colorrange = (0, 1000)) -hmS = heatmap!(axS, stats.Jˢ.max, colormap = :magma, colorrange = (1e-6, 5e-5)) -hmQc = heatmap!(axQc, stats.Qc.max, colormap = :magma, colorrange = (0, 1000)) -hmQv = heatmap!(axQv, stats.Qv.max, colormap = :magma, colorrange = (0, 1000)) -hmu = heatmap!(axu, stats.τx.max * ρ, colormap = :magma, colorrange = (0.1, 2.0)) -hmv = heatmap!(axv, stats.τy.max * ρ, colormap = :magma, colorrange = (0.1, 2.0)) +hmQ = heatmap!(axQ, stats.Jᵀ.max * ρ * cp, colormap = :magma, colorrange = (0, 1000)) +hmS = heatmap!(axS, stats.Jˢ.max, colormap = :magma, colorrange = (1e-6, 5e-5)) +hmQc = heatmap!(axQc, stats.Qc.max, colormap = :magma, colorrange = (0, 1000)) +hmQv = heatmap!(axQv, stats.Qv.max, colormap = :magma, colorrange = (0, 1000)) +hmu = heatmap!(axu, stats.τx.max * ρ, colormap = :magma, colorrange = (0.1, 2.0)) +hmv = heatmap!(axv, stats.τy.max * ρ, colormap = :magma, colorrange = (0.1, 2.0)) Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") Colorbar(fig[2, 0], hmQc, flipaxis=false, label = "W/m²") -Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") +Colorbar(fig[2, 3], hmQv, flipaxis=true, label = "W/m²") Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") -save("fluxes_max.png", fig) +save("fluxes_max.png", fig; px_per_unit) fig = Figure(size = (1200, 900)) -axQ = Axis(fig[1, 1], title = "Heat Flux") -axS = Axis(fig[1, 2], title = "Salinity Flux") + +axQ = Axis(fig[1, 1], title = "Heat Flux") +axS = Axis(fig[1, 2], title = "Salinity Flux") axQc = Axis(fig[2, 1], title = "Sensible Heat Flux") axQv = Axis(fig[2, 2], title = "Latent Heat Flux") -axu = Axis(fig[3, 1], title = "Zonal stress") -axv = Axis(fig[3, 2], title = "Meridional stress") +axu = Axis(fig[3, 1], title = "Zonal stress") +axv = Axis(fig[3, 2], title = "Meridional stress") -hmQ = heatmap!(axQ, stats.Jᵀ.min * ρ * cp, colormap = :haline, colorrange = (-1500, 0)) -hmS = heatmap!(axS, stats.Jˢ.min, colormap = :haline, colorrange = (-5e-6, 0)) -hmQc = heatmap!(axQc, stats.Qc.min, colormap = :haline, colorrange = (-1000, 0)) -hmQv = heatmap!(axQv, stats.Qv.min, colormap = :haline, colorrange = (-1000, 0)) -hmu = heatmap!(axu, stats.τx.min * ρ, colormap = :haline, colorrange = (-2.0, 0.1)) -hmv = heatmap!(axv, stats.τy.min * ρ, colormap = :haline, colorrange = (-2.0, 0.1)) +hmQ = heatmap!(axQ, stats.Jᵀ.min * ρ * cp, colormap = :haline, colorrange = (-1500, 0)) +hmS = heatmap!(axS, stats.Jˢ.min, colormap = :haline, colorrange = (-5e-6, 0)) +hmQc = heatmap!(axQc, stats.Qc.min, colormap = :haline, colorrange = (-1000, 0)) +hmQv = heatmap!(axQv, stats.Qv.min, colormap = :haline, colorrange = (-1000, 0)) +hmu = heatmap!(axu, stats.τx.min * ρ, colormap = :haline, colorrange = (-2.0, 0.1)) +hmv = heatmap!(axv, stats.τy.min * ρ, colormap = :haline, colorrange = (-2.0, 0.1)) Colorbar(fig[1, 0], hmQ, flipaxis=false, label = "W/m²") Colorbar(fig[1, 3], hmS, flipaxis=true, label = "kg/m²/s") @@ -100,4 +106,4 @@ Colorbar(fig[2, 3], hmQv, flipaxis=false, label = "W/m²") Colorbar(fig[3, 0], hmu, flipaxis=false, label = "N/m²") Colorbar(fig[3, 3], hmv, flipaxis=true, label = "N/m²") -save("fluxes_min.png", fig) +save("fluxes_min.png", fig; px_per_unit) From 0d39d05e5393acf9b571c0399dda06dc4e6ad1ae Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 May 2025 17:59:29 +0300 Subject: [PATCH 11/19] compute std in the kernel --- experiments/flux_climatology/flux_climatology.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index c70fc081..a290db8d 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -65,9 +65,6 @@ function update_stats!(stats::FluxStatistics, flux, iteration) grid = flux.grid arch = architecture(grid) launch!(arch, grid, :xy, _update_stats!, stats, flux, iteration) - - parent(stats.std) .= sqrt.(parent(stats.meansq) .- parent(stats.mean).^2) - return nothing end @@ -77,13 +74,15 @@ end inverse_iteration = 1 / (iteration + 1) # use iterative averaging via - # mean = mean * (1 - 1/iteration) + xn / iteration + # mean_n = (x1 + ... + xn) / n -> + # -> mean_{n+1} = mean_n * (1 - 1/(n+1)) * x_{n+1} / (n+1) @inbounds begin f = flux[i, j, 1] stats.mean[i, j, 1] *= 1 - inverse_iteration stats.mean[i, j, 1] += f * inverse_iteration stats.meansq[i, j, 1] *= 1 - inverse_iteration stats.meansq[i, j, 1] += f^2 * inverse_iteration + stats.std[i, j, 1] = sqrt(stats.meansq[i, j, 1] - stats.mean[i, j, 1]^2) stats.max[i, j, 1] = max(stats.max[i, j, 1], f) stats.min[i, j, 1] = min(stats.min[i, j, 1], f) end From 0ec0bd4f70a8045ad973ea91e55bff1cd6e0d123 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 May 2025 12:27:17 +0200 Subject: [PATCH 12/19] add all this --- .../assemble_net_fluxes.jl | 12 +++++++--- .../component_interfaces.jl | 23 ++++++++++++++----- .../InterfaceComputations/radiation.jl | 3 +++ 3 files changed, 29 insertions(+), 9 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 5542794e..ab6d0493 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -103,9 +103,15 @@ end α = stateindex(atmos_ocean_properties.radiation.α, i, j, kᴺ, grid, time) ϵ = stateindex(atmos_ocean_properties.radiation.ϵ, i, j, kᴺ, grid, time) Qu = upwelling_radiation(Tₛ, σ, ϵ) - Qr = (; Qs, Qℓ) - Qd = net_downwelling_radiation(Qr, α, ϵ) - ΣQao = Qd + Qu + Qc + Qv + Qdℓ = downwelling_longwave_radiation(Qℓ, ϵ) + Qds = downwelling_shortwave_radiation(Qs, α) + ΣQao = Qu + Qc + Qv + Qdℓ + Qds + + @inbounds begin + atmos_ocean_fluxes.upwelling_longwave[i, j, 1] = Qu + atmos_ocean_fluxes.downwelling_longwave[i, j, 1] = Qdℓ + atmos_ocean_fluxes.downwelling_shortwave[i, j, 1] = Qds + end # Convert from a mass flux to a volume flux (aka velocity) # by dividing with the density of freshwater. diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl index 73e1203b..35dbe637 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl @@ -178,12 +178,23 @@ function atmosphere_ocean_interface(atmos, velocity_formulation, specific_humidity_formulation) - water_vapor = Field{Center, Center, Nothing}(ocean.model.grid) - latent_heat = Field{Center, Center, Nothing}(ocean.model.grid) - sensible_heat = Field{Center, Center, Nothing}(ocean.model.grid) - x_momentum = Field{Center, Center, Nothing}(ocean.model.grid) - y_momentum = Field{Center, Center, Nothing}(ocean.model.grid) - ao_fluxes = (; latent_heat, sensible_heat, water_vapor, x_momentum, y_momentum) + water_vapor = Field{Center, Center, Nothing}(ocean.model.grid) + latent_heat = Field{Center, Center, Nothing}(ocean.model.grid) + sensible_heat = Field{Center, Center, Nothing}(ocean.model.grid) + x_momentum = Field{Center, Center, Nothing}(ocean.model.grid) + y_momentum = Field{Center, Center, Nothing}(ocean.model.grid) + upwelling_longwave = Field{Center, Center, Nothing}(ocean.model.grid) + downwelling_longwave = Field{Center, Center, Nothing}(ocean.model.grid) + downwelling_shortwave = Field{Center, Center, Nothing}(ocean.model.grid) + + ao_fluxes = (; latent_heat, + sensible_heat, + water_vapor, + x_momentum, + y_momentum, + upwelling_longwave, + downwelling_longwave, + downwelling_shortwave) σ = radiation.stefan_boltzmann_constant αₐₒ = radiation.reflection.ocean diff --git a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl index 88b3ceb5..745306fb 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl @@ -94,3 +94,6 @@ end @inline net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ) = - (1 - α) * Qs - ϵ * Qℓ @inline net_downwelling_radiation(r, α, ϵ) = - (1 - α) * r.Qs - ϵ * r.Qℓ +# Split the individual bands +@inline downwelling_longwave_radiation(Qℓ, ϵ) = - ϵ * Qℓ +@inline downwelling_shortwave_radiation(Qℓ, ϵ) = - (1 - α) * Qs \ No newline at end of file From b5dfad26f68c53863b1980cd864a4e21d312f45e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 May 2025 12:31:19 +0200 Subject: [PATCH 13/19] easier flux computation by storing the flux --- .../flux_climatology/flux_climatology.jl | 29 ++++++++++--------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index a290db8d..0ac4ed4e 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -26,6 +26,7 @@ ClimaOcean.DataWrangling.dataset_defaults.FloatType = Float64 ##### struct FluxStatistics{F} + flux :: F mean :: F meansq :: F std :: F @@ -46,29 +47,31 @@ function FluxStatistics(f::Field) fill!(max, 0) fill!(min, 0) - return FluxStatistics(mean, meansq, std, max, min) + return FluxStatistics(f, mean, meansq, std, max, min) end -Adapt.adapt_structure(to, f::FluxStatistics) = FluxStatistics(Adapt.adapt(to, f.mean), +Adapt.adapt_structure(to, f::FluxStatistics) = FluxStatistics(Adapt.adapt(to, f.flux), + Adapt.adapt(to, f.mean), Adapt.adapt(to, f.meansq), Adapt.adapt(to, f.std), Adapt.adapt(to, f.max), Adapt.adapt(to, f.min)) -on_architecture(arch, f::FluxStatistics) = FluxStatistics(on_architecture(arch, f.mean), +on_architecture(arch, f::FluxStatistics) = FluxStatistics(on_architecture(arch, f.flux), + on_architecture(arch, f.mean), on_architecture(arch, f.meansq), on_architecture(arch, f.std), on_architecture(arch, f.max), on_architecture(arch, f.min)) -function update_stats!(stats::FluxStatistics, flux, iteration) +function update_stats!(stats::FluxStatistics, iteration) grid = flux.grid arch = architecture(grid) - launch!(arch, grid, :xy, _update_stats!, stats, flux, iteration) + launch!(arch, grid, :xy, _update_stats!, stats, iteration) return nothing end -@kernel function _update_stats!(stats, flux, iteration) +@kernel function _update_stats!(stats, iteration) i, j = @index(Global, NTuple) inverse_iteration = 1 / (iteration + 1) @@ -77,7 +80,7 @@ end # mean_n = (x1 + ... + xn) / n -> # -> mean_{n+1} = mean_n * (1 - 1/(n+1)) * x_{n+1} / (n+1) @inbounds begin - f = flux[i, j, 1] + f = stats.flux[i, j, 1] stats.mean[i, j, 1] *= 1 - inverse_iteration stats.mean[i, j, 1] += f * inverse_iteration stats.meansq[i, j, 1] *= 1 - inverse_iteration @@ -107,12 +110,12 @@ function compute_flux_climatology(earth) function update_flux_stats!(earth) iteration = earth.model.clock.iteration - update_stats!(τx, net_fluxes.u, iteration) - update_stats!(τy, net_fluxes.v, iteration) - update_stats!(Jᵀ, net_fluxes.T, iteration) - update_stats!(Jˢ, net_fluxes.S, iteration) - update_stats!(Qc, atmos_ocean_fluxes.sensible_heat, iteration) - update_stats!(Qv, atmos_ocean_fluxes.latent_heat, iteration) + update_stats!(τx, iteration) + update_stats!(τy, iteration) + update_stats!(Jᵀ, iteration) + update_stats!(Jˢ, iteration) + update_stats!(Qc, iteration) + update_stats!(Qv, iteration) return nothing end From 1e048c9162ae51d6b892aac19dbda7c96f6e37e5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 May 2025 12:34:07 +0200 Subject: [PATCH 14/19] improve stats --- .../flux_climatology/flux_climatology.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 0ac4ed4e..f9e40c19 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -65,12 +65,19 @@ on_architecture(arch, f::FluxStatistics) = FluxStatistics(on_architecture(arch, on_architecture(arch, f.min)) function update_stats!(stats::FluxStatistics, iteration) - grid = flux.grid + grid = stats.flux.grid arch = architecture(grid) launch!(arch, grid, :xy, _update_stats!, stats, iteration) return nothing end +function update_stats!(stats::NamedTuple, iteration) + for stat in stats + update_stats!(stat, iteration) + end + return nothing +end + @kernel function _update_stats!(stats, iteration) i, j = @index(Global, NTuple) @@ -110,13 +117,7 @@ function compute_flux_climatology(earth) function update_flux_stats!(earth) iteration = earth.model.clock.iteration - update_stats!(τx, iteration) - update_stats!(τy, iteration) - update_stats!(Jᵀ, iteration) - update_stats!(Jˢ, iteration) - update_stats!(Qc, iteration) - update_stats!(Qv, iteration) - + update_stats!(stats, iteration) return nothing end From 455a12a449f7b3fd46e460180717fa954a78fcc7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 May 2025 14:07:52 +0200 Subject: [PATCH 15/19] fix the edson stability function --- .../similarity_theory_turbulent_fluxes.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl index c2f4e9f0..50fdd814 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl @@ -90,7 +90,7 @@ Keyword Arguments function SimilarityTheoryFluxes(FT::DataType = Oceananigans.defaults.FloatType; von_karman_constant = 0.4, turbulent_prandtl_number = 1, - gustiness_parameter = 1, + gustiness_parameter = 6.5, stability_functions = edson_stability_functions(FT), roughness_lengths = default_roughness_lengths(FT), similarity_form = LogarithmicSimilarityProfile(), @@ -302,7 +302,7 @@ stability function for _stable_ or _unstable_ atmospheric conditions, respective Aˢ :: FT = 0.35 Bˢ :: FT = 0.7 Cˢ :: FT = 0.75 - Dˢ :: FT = 5/0.35 + Dˢ :: FT = 5 / 0.35 Aᵘ :: FT = 15.0 Bᵘ :: FT = 2.0 Cᵘ :: FT = π/2 @@ -326,10 +326,10 @@ end ζ⁻ = min(zero(ζ), ζ) ζ⁺ = max(zero(ζ), ζ) - dζ = min(ζmax, Aˢ * ζ⁺) + dζ = min(ζmax, Aˢ * ζ⁺) # ok # Stability parameter for _stable_ atmospheric conditions - ψₛ = - (Bˢ * ζ⁺ + Cˢ * (ζ⁺ - Dˢ)) * exp(- dζ) - Cˢ * Dˢ + ψₛ = - Bˢ * ζ⁺ - Cˢ * (ζ⁺ - Dˢ) * exp(- dζ) - Cˢ * Dˢ # Stability parameter for _unstable_ atmospheric conditions fᵤ₁ = sqrt(sqrt(1 - Aᵘ * ζ⁻)) From 60969bc6cb9238cad02d5ce44db65bb85d696cbe Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 May 2025 14:11:34 +0200 Subject: [PATCH 16/19] some fixes --- src/DataWrangling/ECCO/ECCO_metadata.jl | 4 ++-- src/OceanSeaIceModels/InterfaceComputations/radiation.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DataWrangling/ECCO/ECCO_metadata.jl b/src/DataWrangling/ECCO/ECCO_metadata.jl index f977004d..4bd45175 100644 --- a/src/DataWrangling/ECCO/ECCO_metadata.jl +++ b/src/DataWrangling/ECCO/ECCO_metadata.jl @@ -127,7 +127,7 @@ ECCO4_short_names = Dict( :net_heat_flux => "oceQnet", :sensible_heat_flux => "EXFhs", :latent_heat_flux => "EXFhl", - :upwelling_longwave => "EXFlwnet" + :net_longwave => "EXFlwnet", :downwelling_shortwave => "oceQsw", :downwelling_longwave => "EXFlwdn", ) @@ -152,7 +152,7 @@ ECCO_location = Dict( :net_heat_flux => (Center, Center, Nothing), :sensible_heat_flux => (Center, Center, Nothing), :latent_heat_flux => (Center, Center, Nothing), - :upwelling_longwave => (Center, Center, Nothing), + :net_longwave => (Center, Center, Nothing), :downwelling_longwave => (Center, Center, Nothing), :downwelling_shortwave => (Center, Center, Nothing), :u_velocity => (Face, Center, Center), diff --git a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl index 745306fb..195af919 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/radiation.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/radiation.jl @@ -96,4 +96,4 @@ end # Split the individual bands @inline downwelling_longwave_radiation(Qℓ, ϵ) = - ϵ * Qℓ -@inline downwelling_shortwave_radiation(Qℓ, ϵ) = - (1 - α) * Qs \ No newline at end of file +@inline downwelling_shortwave_radiation(Qs, α) = - (1 - α) * Qs \ No newline at end of file From ac25d0a1ab7339f24595bed39fc506bece276fa4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 May 2025 14:15:10 +0200 Subject: [PATCH 17/19] add some fluxes --- .../flux_climatology/flux_climatology.jl | 47 ++++++++++++++++--- 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index f9e40c19..0ac50cf1 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -88,11 +88,10 @@ end # -> mean_{n+1} = mean_n * (1 - 1/(n+1)) * x_{n+1} / (n+1) @inbounds begin f = stats.flux[i, j, 1] - stats.mean[i, j, 1] *= 1 - inverse_iteration - stats.mean[i, j, 1] += f * inverse_iteration + stats.mean[i, j, 1] *= 1 - inverse_iteration + stats.mean[i, j, 1] += f * inverse_iteration stats.meansq[i, j, 1] *= 1 - inverse_iteration stats.meansq[i, j, 1] += f^2 * inverse_iteration - stats.std[i, j, 1] = sqrt(stats.meansq[i, j, 1] - stats.mean[i, j, 1]^2) stats.max[i, j, 1] = max(stats.max[i, j, 1], f) stats.min[i, j, 1] = min(stats.min[i, j, 1], f) end @@ -112,8 +111,11 @@ function compute_flux_climatology(earth) atmos_ocean_fluxes = earth.model.interfaces.atmosphere_ocean_interface.fluxes Qc = FluxStatistics(atmos_ocean_fluxes.sensible_heat) Qv = FluxStatistics(atmos_ocean_fluxes.latent_heat) + Qu = FluxStatistics(atmos_ocean_fluxes.upwelling_longwave) + Qs = FluxStatistics(atmos_ocean_fluxes.downwelling_shortwave) + Qℓ = FluxStatistics(atmos_ocean_fluxes.downwelling_longwave) - stats = (; τx, τy, Jᵀ, Jˢ, Qc, Qv) + stats = (; τx, τy, Jᵀ, Jˢ, Qc, Qv, Qu, Qs, Qℓ) function update_flux_stats!(earth) iteration = earth.model.clock.iteration @@ -163,9 +165,9 @@ dataset = ECCO4Monthly() arch = CPU() start_date = DateTime(1992, 1, 1) -end_date = DateTime(1992, 12, 31) +end_date = DateTime(1992, 1, 2) -stop_time = Day(end_date - start_date).value * Oceananigans.Units.days +stop_time = Day(end_date - start_date).value * Oceananigans.Units.days time_indices_in_memory = 13 @@ -179,11 +181,16 @@ v = FieldTimeSeries(v_meta, arch; time_indices_in_memory) T = FieldTimeSeries(T_meta, arch; time_indices_in_memory) S = FieldTimeSeries(S_meta, arch; time_indices_in_memory) -grid = ECCO.ECCO_immersed_grid(arch) +grid = u.grid ocean_model = PrescribedOcean((; u, v, T, S); grid) ocean = Simulation(ocean_model; Δt=3hours, stop_time) +set!(ocean_model.tracers.T, first(T_meta)) +set!(ocean_model.tracers.S, first(S_meta)) +set!(ocean_model.velocities.u, first(u_meta)) +set!(ocean_model.velocities.v, first(v_meta)) + ##### ##### Need to extend a couple of methods ##### @@ -262,3 +269,29 @@ end add_callback!(earth, progress, IterationInterval(16)) stats = compute_flux_climatology(earth) + +Qtecco = Metadata(:net_heat_flux; start_date, end_date, dataset) +Qcecco = Metadata(:sensible_heat_flux; start_date, end_date, dataset) +Qvecco = Metadata(:latent_heat_flux; start_date, end_date, dataset) +Qℓecco = Metadata(:upwelling_longwave; start_date, end_date, dataset) +Qsecco = Metadata(:downwelling_shortwave; start_date, end_date, dataset) + +Qℓecco = FieldTimeSeries(Qℓecco, arch; time_indices_in_memory) +Qtecco = FieldTimeSeries(Qtecco, arch; time_indices_in_memory) +Qcecco = FieldTimeSeries(Qcecco, arch; time_indices_in_memory) +Qvecco = FieldTimeSeries(Qvecco, arch; time_indices_in_memory) +Qsecco = FieldTimeSeries(Qsecco, arch; time_indices_in_memory) + +Q̄tecco = Field{Center, Center, Nothing}(Qtecco.grid) +Q̄cecco = Field{Center, Center, Nothing}(Qtecco.grid) +Q̄vecco = Field{Center, Center, Nothing}(Qtecco.grid) +Q̄ℓecco = Field{Center, Center, Nothing}(Qtecco.grid) +Q̄secco = Field{Center, Center, Nothing}(Qtecco.grid) + +for i in 1:length(Qtecco.times) + Q̄tecco .+= Qtecco[i] ./ length(Qtecco.times) + Q̄cecco .+= Qcecco[i] ./ length(Qcecco.times) + Q̄vecco .+= Qvecco[i] ./ length(Qvecco.times) + Q̄ℓecco .+= Qℓecco[i] ./ length(Qℓecco.times) + Q̄secco .+= Qsecco[i] ./ length(Qsecco.times) +end \ No newline at end of file From d47961f4440c577afc8b76a3df2463c06d9d1fb4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 8 May 2025 15:51:39 +0200 Subject: [PATCH 18/19] small change to trigger docs build --- experiments/flux_climatology/flux_climatology.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 0ac50cf1..ee6295ca 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -109,7 +109,7 @@ function compute_flux_climatology(earth) Jˢ = FluxStatistics(net_fluxes.S) atmos_ocean_fluxes = earth.model.interfaces.atmosphere_ocean_interface.fluxes - Qc = FluxStatistics(atmos_ocean_fluxes.sensible_heat) + Qc = FluxStatistics(atmos_ocean_fluxes.sensible_heat) Qv = FluxStatistics(atmos_ocean_fluxes.latent_heat) Qu = FluxStatistics(atmos_ocean_fluxes.upwelling_longwave) Qs = FluxStatistics(atmos_ocean_fluxes.downwelling_shortwave) From 50861af394c87a31af709ef0102390101cde91b5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 8 May 2025 15:54:16 +0200 Subject: [PATCH 19/19] Update interpolate_atmospheric_state.jl --- .../InterfaceComputations/interpolate_atmospheric_state.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl index 28d0d5f4..6997df68 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl @@ -111,8 +111,8 @@ function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosph # # TODO: find a better design for this that doesn't have redundant # arrays for the barotropic potential - u_potential = forcing_barotropic_potential(ocean) - v_potential = forcing_barotropic_potential(ocean) + u_potential = forcing_barotropic_potential(ocean.model.forcing.u) + v_potential = forcing_barotropic_potential(ocean.model.forcing.v) ρₒ = coupled_model.interfaces.ocean_properties.reference_density if !isnothing(u_potential)