From 5dcbf1df8a6e0c14cb593f017652595588450980 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:06:47 -0600 Subject: [PATCH 01/14] Add a prescribed SST example --- examples/prescribed_sst.jl | 202 ++++++++++++++++++++++++++++++ src/atmospheric_thermodynamics.jl | 8 ++ 2 files changed, 210 insertions(+) create mode 100644 examples/prescribed_sst.jl diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl new file mode 100644 index 00000000..7a992616 --- /dev/null +++ b/examples/prescribed_sst.jl @@ -0,0 +1,202 @@ +using Oceananigans +using Oceananigans.Units +using Printf +using AquaSkyLES + +Nx = Nz = 128 +Ny = 1 +Lz = 4 * 1024 +grid = RectilinearGrid(size = (Nx, Ny, Nz), + x = (0, 2Lz), + y = (0, 2Lz), + z = (0, Lz), + topology = (Periodic, Periodic, Bounded)) + +p₀ = 101325 # Pa +θ₀ = 288 # K +reference_state = AquaSkyLES.ReferenceState(base_pressure=p₀, potential_temperature=θ₀) +buoyancy = AquaSkyLES.MoistAirBuoyancy(; reference_state) +thermodynamics = buoyancy.thermodynamics +condensation = thermodynamics.condensation + +ρ₀ = AquaSkyLES.base_density(buoyancy) # air density at z=0 +cₚ = buoyancy.thermodynamics.dry_air.heat_capacity + +parameters = (; + drag_coefficient = 2e-3, + heat_transfer_coefficient = 1e-3, + vapor_transfer_coefficient = 1e-3, + sea_surface_temperature = θ₀ + 10, + reference_state, + thermodynamics, + condensation +) + +@inline surface_saturation_specific_humidity(T, parameters) = + AquaSkyLES.saturation_specific_humidity(T, zero(T), + parameters.reference_state, + parameters.thermodynamics, + parameters.condensation) + +@inline function friction_velocity(x, y, t, u, v, parameters) + Cᴰ = parameters.drag_coefficient + return sqrt(Cᴰ * (u^2 + v^2)) +end + +@inline function x_momentum_flux(x, y, t, u, v, parameters) + u★ = friction_velocity(x, y, t, u, v, parameters) + return - u★^2 * u / sqrt(u^2 + v^2) +end + +@inline function y_momentum_flux(x, y, t, u, v, parameters) + u★ = friction_velocity(x, y, t, u, v, parameters) + return - u★^2 * v / sqrt(u^2 + v^2) +end + +@inline function temperature_flux(x, y, t, u, v, θ, parameters) + u★ = friction_velocity(x, y, t, u, v, parameters) + θˢ = parameters.sea_surface_temperature + Cᴰ = parameters.drag_coefficient + Cᴴ = parameters.heat_transfer_coefficient + Δθ = θ - θˢ + θ★ = Cᴴ / sqrt(Cᴰ) * Δθ + return u★ * θ★ +end + +@inline function vapor_flux(x, y, t, u, v, q, parameters) + u★ = friction_velocity(x, y, t, u, v, parameters) + # Note: typically we would compute this using the saturation specific humidity + # given the sea surface temperature + # θˢ = parameters.sea_surface_temperature + # qˢ = surface_saturation_specific_humidity(θˢ, parameters) + qˢ = zero(u★) + Cᴰ = parameters.drag_coefficient + Cᵛ = parameters.vapor_transfer_coefficient + Δq = q - qˢ + q★ = Cᵛ / sqrt(Cᴰ) * Δq + return u★ * q★ +end + +u_surface_flux = FluxBoundaryCondition(x_momentum_flux; field_dependencies=(:u, :v), parameters) +v_surface_flux = FluxBoundaryCondition(y_momentum_flux; field_dependencies=(:u, :v), parameters) +θ_surface_flux = FluxBoundaryCondition(temperature_flux; field_dependencies=(:u, :v, :θ), parameters) +q_surface_flux = FluxBoundaryCondition(vapor_flux; field_dependencies=(:u, :v, :q), parameters) + +u_bcs = FieldBoundaryConditions(bottom=u_surface_flux) +v_bcs = FieldBoundaryConditions(bottom=v_surface_flux) +θ_bcs = FieldBoundaryConditions(bottom=θ_surface_flux) +q_bcs = FieldBoundaryConditions(bottom=q_surface_flux) + +advection = WENO() #(momentum=WENO(), θ=WENO(), q=WENO(bounds=(0, 1))) +tracers = (:θ, :q) +model = NonhydrostaticModel(; grid, advection, buoyancy, + tracers = (:θ, :q), + boundary_conditions = (u=u_bcs, v=v_bcs, θ=θ_bcs, q=q_bcs)) + +Lz = grid.Lz +Δθ = 5 # K +Tₛ = reference_state.θ # K +θᵢ(x, y, z) = Tₛ + Δθ * z / Lz + 1e-2 * Δθ * randn() +qᵢ(x, y, z) = 1e-2 + 1e-5 * rand() +set!(model, θ=θᵢ, q=qᵢ) + +simulation = Simulation(model, Δt=10, stop_time=2hours) +conjure_time_step_wizard!(simulation, cfl=0.7) + +T = AquaSkyLES.TemperatureField(model) +qˡ = AquaSkyLES.CondensateField(model, T) +qᵛ★ = AquaSkyLES.SaturationField(model, T) +δ = Field(model.tracers.q - qᵛ★) + +function progress(sim) + compute!(T) + compute!(qˡ) + compute!(δ) + q = sim.model.tracers.q + θ = sim.model.tracers.θ + u, v, w = sim.model.velocities + + umax = maximum(u) + vmax = maximum(v) + wmax = maximum(w) + + qmin = minimum(q) + qmax = maximum(q) + qˡmax = maximum(qˡ) + δmax = maximum(δ) + + θmin = minimum(θ) + θmax = maximum(θ) + + msg = @sprintf("Iter: %d, t = %s, max|u|: (%.2e, %.2e, %.2e)", + iteration(sim), prettytime(sim), umax, vmax, wmax) + + msg *= @sprintf(", extrema(q): (%.2e, %.2e), max(qˡ): %.2e, min(δ): %.2e, extrema(θ): (%.2e, %.2e)", + qmin, qmax, qˡmax, δmax, θmin, θmax) + + @info msg + + return nothing +end + +add_callback!(simulation, progress, IterationInterval(10)) + +outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) + +ow = JLD2Writer(model, outputs, + filename = "free_convection.jld2", + schedule = IterationInterval(10), + overwrite_existing = true) + +simulation.output_writers[:jld2] = ow + +run!(simulation) + +θt = FieldTimeSeries("free_convection.jld2", "θ") +Tt = FieldTimeSeries("free_convection.jld2", "T") +qt = FieldTimeSeries("free_convection.jld2", "q") +qˡt = FieldTimeSeries("free_convection.jld2", "qˡ") +times = qt.times +Nt = length(θt) + +using GLMakie, Printf + +n = Observable(length(θt)) + +θn = @lift θt[$n] +qn = @lift qt[$n] +Tn = @lift Tt[$n] +qˡn = @lift qˡt[$n] +title = @lift "t = $(prettytime(times[$n]))" + +fig = Figure(size=(800, 400), fontsize=12) +axθ = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)") +axq = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)") +axT = Axis(fig[2, 1], xlabel="x (m)", ylabel="z (m)") +axqˡ = Axis(fig[2, 2], xlabel="x (m)", ylabel="z (m)") + +fig[0, :] = Label(fig, title, fontsize=22, tellwidth=false) + +Tmin = minimum(Tt) +Tmax = maximum(Tt) + +hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ)) +hmq = heatmap!(axq, qn, colorrange=(0, 2e-2), colormap=:magma) +hmT = heatmap!(axT, Tn, colorrange=(Tmin, Tmax)) +hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, 2e-4), colormap=:magma) + +# Label(fig[0, 1], "θ", tellwidth=false) +# Label(fig[0, 2], "q", tellwidth=false) +# Label(fig[0, 1], "θ", tellwidth=false) +# Label(fig[0, 2], "q", tellwidth=false) + +Colorbar(fig[1, 0], hmθ, label = "θ [K]", vertical=true) +Colorbar(fig[1, 3], hmq, label = "q", vertical=true) +Colorbar(fig[2, 0], hmT, label = "T [K]", vertical=true) +Colorbar(fig[2, 3], hmqˡ, label = "qˡ", vertical=true) + +fig + +record(fig, "free_convection.mp4", 1:Nt, framerate=12) do nn + n[] = nn +end diff --git a/src/atmospheric_thermodynamics.jl b/src/atmospheric_thermodynamics.jl index 87f093f8..a5607579 100644 --- a/src/atmospheric_thermodynamics.jl +++ b/src/atmospheric_thermodynamics.jl @@ -284,6 +284,14 @@ end return ref.p₀ * (1 - g * z / (cᵖᵈ * ref.θ))^ϰᵈ⁻¹ end +""" + saturation_specific_humidity(T, z, reference_state, thermodynamics, phase_transition) + +Using the Clasius-Claperyon relation for saturation vapor pressure, compute theo +saturation specific humidity at temperature `T`, height `z`, given the +`reference_state`, `thermodynamics` constants, and `phase_transition` constants representing +parameters representing either condensation (vapor to liquid) or deposition (vapor to solid). +""" @inline function saturation_specific_humidity(T, z, ref::ReferenceState, thermo, phase_transition) ρ = reference_density(z, ref, thermo) return saturation_specific_humidity(T, ρ, thermo, phase_transition) From 1d8833827a870c51fda00a12d7716ab352342dcf Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:31:49 -0600 Subject: [PATCH 02/14] handle U=0 --- examples/prescribed_sst.jl | 51 ++++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 7a992616..9cde8865 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -3,7 +3,8 @@ using Oceananigans.Units using Printf using AquaSkyLES -Nx = Nz = 128 +Nx = 4 +Nz = 128 Ny = 1 Lz = 4 * 1024 grid = RectilinearGrid(size = (Nx, Ny, Nz), @@ -19,38 +20,45 @@ buoyancy = AquaSkyLES.MoistAirBuoyancy(; reference_state) thermodynamics = buoyancy.thermodynamics condensation = thermodynamics.condensation -ρ₀ = AquaSkyLES.base_density(buoyancy) # air density at z=0 -cₚ = buoyancy.thermodynamics.dry_air.heat_capacity +# ρ₀ = AquaSkyLES.base_density(buoyancy) # air density at z=0 +# cₚ = buoyancy.thermodynamics.dry_air.heat_capacity + +Δz = Lz / Nz # grid spacing parameters = (; - drag_coefficient = 2e-3, - heat_transfer_coefficient = 1e-3, - vapor_transfer_coefficient = 1e-3, - sea_surface_temperature = θ₀ + 10, - reference_state, + drag_coefficient = 1e-4, + heat_transfer_coefficient = 1e-5, + vapor_transfer_coefficient = 1e-5, + sea_surface_temperature = θ₀ + 1, + base_air_density = AquaSkyLES.base_density(buoyancy), # air density at z=0, thermodynamics, condensation ) +# Utility for computing the saturation specific humidity at the sea surface @inline surface_saturation_specific_humidity(T, parameters) = - AquaSkyLES.saturation_specific_humidity(T, zero(T), - parameters.reference_state, + AquaSkyLES.saturation_specific_humidity(T, parameters.base_air_density, parameters.thermodynamics, parameters.condensation) @inline function friction_velocity(x, y, t, u, v, parameters) Cᴰ = parameters.drag_coefficient - return sqrt(Cᴰ * (u^2 + v^2)) + Δu = u # stationary ocean + Δv = v # stationary ocean + return sqrt(Cᴰ * (Δu^2 + Δv^2)) end +# Take care to handle U = 0 @inline function x_momentum_flux(x, y, t, u, v, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - return - u★^2 * u / sqrt(u^2 + v^2) + U = sqrt(u^2 + v^2) + return - u★^2 * u / U * (U > 0) end @inline function y_momentum_flux(x, y, t, u, v, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - return - u★^2 * v / sqrt(u^2 + v^2) + U = sqrt(u^2 + v^2) + return - u★^2 * v / U * (U > 0) end @inline function temperature_flux(x, y, t, u, v, θ, parameters) @@ -59,22 +67,21 @@ end Cᴰ = parameters.drag_coefficient Cᴴ = parameters.heat_transfer_coefficient Δθ = θ - θˢ + # Using the scaling argument: u★ θ★ = Cᴴ * U * Δθ θ★ = Cᴴ / sqrt(Cᴰ) * Δθ - return u★ * θ★ + return - u★ * θ★ end @inline function vapor_flux(x, y, t, u, v, q, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - # Note: typically we would compute this using the saturation specific humidity - # given the sea surface temperature - # θˢ = parameters.sea_surface_temperature - # qˢ = surface_saturation_specific_humidity(θˢ, parameters) - qˢ = zero(u★) + θˢ = parameters.sea_surface_temperature + qˢ = surface_saturation_specific_humidity(θˢ, parameters) Cᴰ = parameters.drag_coefficient Cᵛ = parameters.vapor_transfer_coefficient Δq = q - qˢ + # Using the scaling argument: u★ q★ = Cᵛ * U * Δq q★ = Cᵛ / sqrt(Cᴰ) * Δq - return u★ * q★ + return - u★ * q★ end u_surface_flux = FluxBoundaryCondition(x_momentum_flux; field_dependencies=(:u, :v), parameters) @@ -139,12 +146,12 @@ function progress(sim) return nothing end -add_callback!(simulation, progress, IterationInterval(10)) +add_callback!(simulation, progress, IterationInterval(1)) outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) ow = JLD2Writer(model, outputs, - filename = "free_convection.jld2", + filename = "prescribed_sst_convection.jld2", schedule = IterationInterval(10), overwrite_existing = true) From abdb6889f76eac3724600603eeb4debda9247557 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:35:19 -0600 Subject: [PATCH 03/14] change movie name and boost resolution --- examples/prescribed_sst.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 9cde8865..297199e7 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -3,7 +3,7 @@ using Oceananigans.Units using Printf using AquaSkyLES -Nx = 4 +Nx = 128 Nz = 128 Ny = 1 Lz = 4 * 1024 @@ -204,6 +204,6 @@ Colorbar(fig[2, 3], hmqˡ, label = "qˡ", vertical=true) fig -record(fig, "free_convection.mp4", 1:Nt, framerate=12) do nn +record(fig, "prescribed_sst.mp4", 1:Nt, framerate=12) do nn n[] = nn end From 9089b00329a229601974b01d9883f167fbaad73b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:43:56 -0600 Subject: [PATCH 04/14] add gustiness --- examples/prescribed_sst.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 297199e7..0228b9ef 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -30,6 +30,7 @@ parameters = (; heat_transfer_coefficient = 1e-5, vapor_transfer_coefficient = 1e-5, sea_surface_temperature = θ₀ + 1, + gust_speed = 1e-2, base_air_density = AquaSkyLES.base_density(buoyancy), # air density at z=0, thermodynamics, condensation @@ -45,19 +46,19 @@ parameters = (; Cᴰ = parameters.drag_coefficient Δu = u # stationary ocean Δv = v # stationary ocean - return sqrt(Cᴰ * (Δu^2 + Δv^2)) + return sqrt(Cᴰ * (Δu^2 + Δv^2)) + 1e-2 end # Take care to handle U = 0 @inline function x_momentum_flux(x, y, t, u, v, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - U = sqrt(u^2 + v^2) + U = sqrt(u^2 + v^2) + parameters.gust_speed return - u★^2 * u / U * (U > 0) end @inline function y_momentum_flux(x, y, t, u, v, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - U = sqrt(u^2 + v^2) + U = sqrt(u^2 + v^2) + parameters.gust_speed return - u★^2 * v / U * (U > 0) end @@ -80,7 +81,7 @@ end Cᵛ = parameters.vapor_transfer_coefficient Δq = q - qˢ # Using the scaling argument: u★ q★ = Cᵛ * U * Δq - q★ = Cᵛ / sqrt(Cᴰ) * Δq + q★ = Cᵛ / sqrt(Cᴰ) * Δq return - u★ * q★ end From dcb8be003fa8e06d3dd61b28958d2764af9cca32 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:45:07 -0600 Subject: [PATCH 05/14] less frequent progress --- examples/prescribed_sst.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 0228b9ef..31100543 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -147,7 +147,7 @@ function progress(sim) return nothing end -add_callback!(simulation, progress, IterationInterval(1)) +add_callback!(simulation, progress, IterationInterval(10)) outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) From 51e3d1c7a4b85bb1a0a00b90c7d8aae8d37eb988 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:48:30 -0600 Subject: [PATCH 06/14] bigger coefficients --- examples/prescribed_sst.jl | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 31100543..d1abc3f5 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -26,11 +26,11 @@ condensation = thermodynamics.condensation Δz = Lz / Nz # grid spacing parameters = (; - drag_coefficient = 1e-4, - heat_transfer_coefficient = 1e-5, - vapor_transfer_coefficient = 1e-5, - sea_surface_temperature = θ₀ + 1, - gust_speed = 1e-2, + drag_coefficient = 1e-3, + heat_transfer_coefficient = 1e-3, + vapor_transfer_coefficient = 5e-4, + sea_surface_temperature = θ₀ + 10, + gust_speed = 1e-1, base_air_density = AquaSkyLES.base_density(buoyancy), # air density at z=0, thermodynamics, condensation @@ -46,19 +46,19 @@ parameters = (; Cᴰ = parameters.drag_coefficient Δu = u # stationary ocean Δv = v # stationary ocean - return sqrt(Cᴰ * (Δu^2 + Δv^2)) + 1e-2 + return sqrt(Cᴰ * (Δu^2 + Δv^2)) + parameters.gust_speed end # Take care to handle U = 0 @inline function x_momentum_flux(x, y, t, u, v, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - U = sqrt(u^2 + v^2) + parameters.gust_speed + U = sqrt(u^2 + v^2) return - u★^2 * u / U * (U > 0) end @inline function y_momentum_flux(x, y, t, u, v, parameters) u★ = friction_velocity(x, y, t, u, v, parameters) - U = sqrt(u^2 + v^2) + parameters.gust_speed + U = sqrt(u^2 + v^2) return - u★^2 * v / U * (U > 0) end @@ -70,7 +70,8 @@ end Δθ = θ - θˢ # Using the scaling argument: u★ θ★ = Cᴴ * U * Δθ θ★ = Cᴴ / sqrt(Cᴰ) * Δθ - return - u★ * θ★ + @show Jθ = - u★ * θ★ + return Jθ end @inline function vapor_flux(x, y, t, u, v, q, parameters) @@ -82,7 +83,8 @@ end Δq = q - qˢ # Using the scaling argument: u★ q★ = Cᵛ * U * Δq q★ = Cᵛ / sqrt(Cᴰ) * Δq - return - u★ * q★ + @show Jq = - u★ * q★ + return Jq end u_surface_flux = FluxBoundaryCondition(x_momentum_flux; field_dependencies=(:u, :v), parameters) From 6795a40252c04a1c9a4b4aa0edf80b370716fdcb Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:53:07 -0600 Subject: [PATCH 07/14] fix typo --- src/atmospheric_thermodynamics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atmospheric_thermodynamics.jl b/src/atmospheric_thermodynamics.jl index a5607579..3ac6fede 100644 --- a/src/atmospheric_thermodynamics.jl +++ b/src/atmospheric_thermodynamics.jl @@ -287,7 +287,7 @@ end """ saturation_specific_humidity(T, z, reference_state, thermodynamics, phase_transition) -Using the Clasius-Claperyon relation for saturation vapor pressure, compute theo +Using the Clasius-Claperyon relation for saturation vapor pressure, compute the saturation specific humidity at temperature `T`, height `z`, given the `reference_state`, `thermodynamics` constants, and `phase_transition` constants representing parameters representing either condensation (vapor to liquid) or deposition (vapor to solid). From 325b5210fe5d3d9c68b2d4093981bcf566c24c1e Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 12 May 2025 18:54:39 -0600 Subject: [PATCH 08/14] run longer --- examples/prescribed_sst.jl | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index d1abc3f5..700f1df6 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -14,7 +14,7 @@ grid = RectilinearGrid(size = (Nx, Ny, Nz), topology = (Periodic, Periodic, Bounded)) p₀ = 101325 # Pa -θ₀ = 288 # K +θ₀ = 285 # K reference_state = AquaSkyLES.ReferenceState(base_pressure=p₀, potential_temperature=θ₀) buoyancy = AquaSkyLES.MoistAirBuoyancy(; reference_state) thermodynamics = buoyancy.thermodynamics @@ -28,9 +28,9 @@ condensation = thermodynamics.condensation parameters = (; drag_coefficient = 1e-3, heat_transfer_coefficient = 1e-3, - vapor_transfer_coefficient = 5e-4, + vapor_transfer_coefficient = 1e-3, sea_surface_temperature = θ₀ + 10, - gust_speed = 1e-1, + gust_speed = 1e-2, base_air_density = AquaSkyLES.base_density(buoyancy), # air density at z=0, thermodynamics, condensation @@ -70,8 +70,7 @@ end Δθ = θ - θˢ # Using the scaling argument: u★ θ★ = Cᴴ * U * Δθ θ★ = Cᴴ / sqrt(Cᴰ) * Δθ - @show Jθ = - u★ * θ★ - return Jθ + return - u★ * θ★ end @inline function vapor_flux(x, y, t, u, v, q, parameters) @@ -83,8 +82,7 @@ end Δq = q - qˢ # Using the scaling argument: u★ q★ = Cᵛ * U * Δq q★ = Cᵛ / sqrt(Cᴰ) * Δq - @show Jq = - u★ * q★ - return Jq + return - u★ * q★ end u_surface_flux = FluxBoundaryCondition(x_momentum_flux; field_dependencies=(:u, :v), parameters) @@ -110,7 +108,7 @@ Tₛ = reference_state.θ # K qᵢ(x, y, z) = 1e-2 + 1e-5 * rand() set!(model, θ=θᵢ, q=qᵢ) -simulation = Simulation(model, Δt=10, stop_time=2hours) +simulation = Simulation(model, Δt=10, stop_time=4hours) conjure_time_step_wizard!(simulation, cfl=0.7) T = AquaSkyLES.TemperatureField(model) @@ -155,17 +153,17 @@ outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) ow = JLD2Writer(model, outputs, filename = "prescribed_sst_convection.jld2", - schedule = IterationInterval(10), + schedule = TimeInterval(2minutes), overwrite_existing = true) simulation.output_writers[:jld2] = ow run!(simulation) -θt = FieldTimeSeries("free_convection.jld2", "θ") -Tt = FieldTimeSeries("free_convection.jld2", "T") -qt = FieldTimeSeries("free_convection.jld2", "q") -qˡt = FieldTimeSeries("free_convection.jld2", "qˡ") +θt = FieldTimeSeries("prescribed_sst_convection.jld2", "θ") +Tt = FieldTimeSeries("prescribed_sst_convection.jld2", "T") +qt = FieldTimeSeries("prescribed_sst_convection.jld2", "q") +qˡt = FieldTimeSeries("prescribed_sst_convection.jld2", "qˡ") times = qt.times Nt = length(θt) From e0c7583d8d0c90add1f0a2e94ea137234d45bace Mon Sep 17 00:00:00 2001 From: Madelaine Gamble Rosevear Date: Tue, 13 May 2025 14:11:20 +1000 Subject: [PATCH 09/14] attempt at outputting fluxes --- examples/prescribed_sst.jl | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 700f1df6..68aaba1a 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -6,7 +6,7 @@ using AquaSkyLES Nx = 128 Nz = 128 Ny = 1 -Lz = 4 * 1024 +Lz = 2 * 1024 grid = RectilinearGrid(size = (Nx, Ny, Nz), x = (0, 2Lz), y = (0, 2Lz), @@ -30,7 +30,7 @@ parameters = (; heat_transfer_coefficient = 1e-3, vapor_transfer_coefficient = 1e-3, sea_surface_temperature = θ₀ + 10, - gust_speed = 1e-2, + gust_speed = 1e-2, # directly added to friction velocity (i.e. not multiplied by drag coefficient Cᴰ) base_air_density = AquaSkyLES.base_density(buoyancy), # air density at z=0, thermodynamics, condensation @@ -102,13 +102,13 @@ model = NonhydrostaticModel(; grid, advection, buoyancy, boundary_conditions = (u=u_bcs, v=v_bcs, θ=θ_bcs, q=q_bcs)) Lz = grid.Lz -Δθ = 5 # K +Δθ = 2.5 # K Tₛ = reference_state.θ # K θᵢ(x, y, z) = Tₛ + Δθ * z / Lz + 1e-2 * Δθ * randn() qᵢ(x, y, z) = 1e-2 + 1e-5 * rand() set!(model, θ=θᵢ, q=qᵢ) -simulation = Simulation(model, Δt=10, stop_time=4hours) +simulation = Simulation(model, Δt=10, stop_time=1hours) conjure_time_step_wizard!(simulation, cfl=0.7) T = AquaSkyLES.TemperatureField(model) @@ -116,6 +116,19 @@ qˡ = AquaSkyLES.CondensateField(model, T) qᵛ★ = AquaSkyLES.SaturationField(model, T) δ = Field(model.tracers.q - qᵛ★) +# new field for outputting Flux BCs +u_flux = Field{Face, Face, Center}(model.grid, indices=(:, :, 1)) + +# something like this? +function update_surface_fluxes!(model, parameters) + u = model.velocities.u[:, :, 1] + v = model.velocities.v[:, :, 1] + + u_flux[:, :] = x_momentum_flux(???, u, v, parameters) + compute!(u_flux) +end + + function progress(sim) compute!(T) compute!(qˡ) @@ -124,6 +137,9 @@ function progress(sim) θ = sim.model.tracers.θ u, v, w = sim.model.velocities + # something like this? + update_surface_fluxes!(sim.model, parameters) + umax = maximum(u) vmax = maximum(v) wmax = maximum(w) @@ -149,7 +165,12 @@ end add_callback!(simulation, progress, IterationInterval(10)) +# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) + outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) +surface_fluxes = (; u_flux) + +outputs = merge(outputs, surface_fluxes) ow = JLD2Writer(model, outputs, filename = "prescribed_sst_convection.jld2", @@ -189,9 +210,9 @@ Tmin = minimum(Tt) Tmax = maximum(Tt) hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ)) -hmq = heatmap!(axq, qn, colorrange=(0, 2e-2), colormap=:magma) +hmq = heatmap!(axq, qn, colorrange=(0.97e-2, 1.05e-2), colormap=:magma) hmT = heatmap!(axT, Tn, colorrange=(Tmin, Tmax)) -hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, 2e-4), colormap=:magma) +hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, 1.5e-3), colormap=:magma) # Label(fig[0, 1], "θ", tellwidth=false) # Label(fig[0, 2], "q", tellwidth=false) From 7eab2b81743ddf1a5657cc949854ec8ed09d30d2 Mon Sep 17 00:00:00 2001 From: Madelaine Gamble Rosevear Date: Wed, 14 May 2025 11:09:56 +1000 Subject: [PATCH 10/14] updated fluxes to discrete form --- examples/prescribed_sst.jl | 63 ++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 36 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 68aaba1a..ebf7af1e 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -42,53 +42,62 @@ parameters = (; parameters.thermodynamics, parameters.condensation) -@inline function friction_velocity(x, y, t, u, v, parameters) + +@inline function friction_velocity(i, j, grid, clock, model_fields, parameters) Cᴰ = parameters.drag_coefficient + u = model_fields.u[i, j, 1] + v = model_fields.v[i, j, 1] Δu = u # stationary ocean Δv = v # stationary ocean return sqrt(Cᴰ * (Δu^2 + Δv^2)) + parameters.gust_speed end # Take care to handle U = 0 -@inline function x_momentum_flux(x, y, t, u, v, parameters) - u★ = friction_velocity(x, y, t, u, v, parameters) +@inline function x_momentum_flux(i, j, grid, clock, model_fields, parameters) + u = model_fields.u[i, j, 1] + v = model_fields.v[i, j, 1] + u★ = friction_velocity(i, j, grid, clock, model_fields, parameters) U = sqrt(u^2 + v^2) return - u★^2 * u / U * (U > 0) end -@inline function y_momentum_flux(x, y, t, u, v, parameters) - u★ = friction_velocity(x, y, t, u, v, parameters) +@inline function y_momentum_flux(i, j, grid, clock, model_fields, parameters) + u = model_fields.u[i, j, 1] + v = model_fields.v[i, j, 1] + u★ = friction_velocity(i, j, grid, clock, model_fields, parameters) U = sqrt(u^2 + v^2) return - u★^2 * v / U * (U > 0) end -@inline function temperature_flux(x, y, t, u, v, θ, parameters) - u★ = friction_velocity(x, y, t, u, v, parameters) +@inline function temperature_flux(i, j, grid, clock, model_fields, parameters) + u★ = friction_velocity(i, j, grid, clock, model_fields, parameters) θˢ = parameters.sea_surface_temperature Cᴰ = parameters.drag_coefficient Cᴴ = parameters.heat_transfer_coefficient - Δθ = θ - θˢ + Δθ = model_fields.θ[i, j, 1] - θˢ # Using the scaling argument: u★ θ★ = Cᴴ * U * Δθ θ★ = Cᴴ / sqrt(Cᴰ) * Δθ return - u★ * θ★ end -@inline function vapor_flux(x, y, t, u, v, q, parameters) - u★ = friction_velocity(x, y, t, u, v, parameters) +@inline function vapor_flux(i, j, grid, clock, model_fields, parameters) + u★ = friction_velocity(i, j, grid, clock, model_fields, parameters) θˢ = parameters.sea_surface_temperature qˢ = surface_saturation_specific_humidity(θˢ, parameters) Cᴰ = parameters.drag_coefficient Cᵛ = parameters.vapor_transfer_coefficient - Δq = q - qˢ + Δq = model_fields.q[i, j, 1] - qˢ # Using the scaling argument: u★ q★ = Cᵛ * U * Δq q★ = Cᵛ / sqrt(Cᴰ) * Δq return - u★ * q★ end -u_surface_flux = FluxBoundaryCondition(x_momentum_flux; field_dependencies=(:u, :v), parameters) -v_surface_flux = FluxBoundaryCondition(y_momentum_flux; field_dependencies=(:u, :v), parameters) -θ_surface_flux = FluxBoundaryCondition(temperature_flux; field_dependencies=(:u, :v, :θ), parameters) -q_surface_flux = FluxBoundaryCondition(vapor_flux; field_dependencies=(:u, :v, :q), parameters) +model_fields = merge(model.velocities, model.tracers) + +u_surface_flux = FluxBoundaryCondition(x_momentum_flux; discrete_form=true, parameters) +v_surface_flux = FluxBoundaryCondition(y_momentum_flux; discrete_form=true, parameters) +θ_surface_flux = FluxBoundaryCondition(temperature_flux; discrete_form=true, parameters) +q_surface_flux = FluxBoundaryCondition(vapor_flux; discrete_form=true, parameters) u_bcs = FieldBoundaryConditions(bottom=u_surface_flux) v_bcs = FieldBoundaryConditions(bottom=v_surface_flux) @@ -108,7 +117,7 @@ Tₛ = reference_state.θ # K qᵢ(x, y, z) = 1e-2 + 1e-5 * rand() set!(model, θ=θᵢ, q=qᵢ) -simulation = Simulation(model, Δt=10, stop_time=1hours) +simulation = Simulation(model, Δt=10, stop_time=4hours) conjure_time_step_wizard!(simulation, cfl=0.7) T = AquaSkyLES.TemperatureField(model) @@ -116,19 +125,6 @@ qˡ = AquaSkyLES.CondensateField(model, T) qᵛ★ = AquaSkyLES.SaturationField(model, T) δ = Field(model.tracers.q - qᵛ★) -# new field for outputting Flux BCs -u_flux = Field{Face, Face, Center}(model.grid, indices=(:, :, 1)) - -# something like this? -function update_surface_fluxes!(model, parameters) - u = model.velocities.u[:, :, 1] - v = model.velocities.v[:, :, 1] - - u_flux[:, :] = x_momentum_flux(???, u, v, parameters) - compute!(u_flux) -end - - function progress(sim) compute!(T) compute!(qˡ) @@ -137,8 +133,8 @@ function progress(sim) θ = sim.model.tracers.θ u, v, w = sim.model.velocities - # something like this? - update_surface_fluxes!(sim.model, parameters) + # # something like this? + # update_surface_fluxes!(sim.model, parameters) umax = maximum(u) vmax = maximum(v) @@ -165,12 +161,7 @@ end add_callback!(simulation, progress, IterationInterval(10)) -# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) - outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) -surface_fluxes = (; u_flux) - -outputs = merge(outputs, surface_fluxes) ow = JLD2Writer(model, outputs, filename = "prescribed_sst_convection.jld2", From 834da39d3f95f6e4e6d7bce7dc6849f7f6ab350e Mon Sep 17 00:00:00 2001 From: Madelaine Gamble Rosevear Date: Wed, 14 May 2025 11:13:10 +1000 Subject: [PATCH 11/14] tidy up --- examples/prescribed_sst.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index ebf7af1e..7b0daa4d 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -133,9 +133,6 @@ function progress(sim) θ = sim.model.tracers.θ u, v, w = sim.model.velocities - # # something like this? - # update_surface_fluxes!(sim.model, parameters) - umax = maximum(u) vmax = maximum(v) wmax = maximum(w) From c2f4da55053628d65b1967c8ba7d466d3e3338df Mon Sep 17 00:00:00 2001 From: Madelaine Gamble Rosevear Date: Wed, 14 May 2025 11:45:16 +1000 Subject: [PATCH 12/14] removed line where I try to define tuple model_fields --- examples/prescribed_sst.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 7b0daa4d..4fc5ec13 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -92,8 +92,6 @@ end return - u★ * q★ end -model_fields = merge(model.velocities, model.tracers) - u_surface_flux = FluxBoundaryCondition(x_momentum_flux; discrete_form=true, parameters) v_surface_flux = FluxBoundaryCondition(y_momentum_flux; discrete_form=true, parameters) θ_surface_flux = FluxBoundaryCondition(temperature_flux; discrete_form=true, parameters) @@ -117,7 +115,7 @@ Tₛ = reference_state.θ # K qᵢ(x, y, z) = 1e-2 + 1e-5 * rand() set!(model, θ=θᵢ, q=qᵢ) -simulation = Simulation(model, Δt=10, stop_time=4hours) +simulation = Simulation(model, Δt=10, stop_time=hours) conjure_time_step_wizard!(simulation, cfl=0.7) T = AquaSkyLES.TemperatureField(model) @@ -125,6 +123,7 @@ qˡ = AquaSkyLES.CondensateField(model, T) qᵛ★ = AquaSkyLES.SaturationField(model, T) δ = Field(model.tracers.q - qᵛ★) + function progress(sim) compute!(T) compute!(qˡ) From bb8f6985d5769a35653b6d4ab94d57eab3417b99 Mon Sep 17 00:00:00 2001 From: Madelaine Gamble Rosevear Date: Wed, 14 May 2025 12:01:15 +1000 Subject: [PATCH 13/14] longer test --- examples/prescribed_sst.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 4fc5ec13..b3f69764 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -115,7 +115,7 @@ Tₛ = reference_state.θ # K qᵢ(x, y, z) = 1e-2 + 1e-5 * rand() set!(model, θ=θᵢ, q=qᵢ) -simulation = Simulation(model, Δt=10, stop_time=hours) +simulation = Simulation(model, Δt=10, stop_time=4hours) conjure_time_step_wizard!(simulation, cfl=0.7) T = AquaSkyLES.TemperatureField(model) From 92f748f87fd4a22eec238c41bc2b14c62014f1c5 Mon Sep 17 00:00:00 2001 From: Madelaine Gamble Rosevear Date: Wed, 28 May 2025 16:36:50 +1000 Subject: [PATCH 14/14] added flux plots --- examples/prescribed_sst.jl | 73 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 70 insertions(+), 3 deletions(-) diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index b3f69764..dc900366 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -123,6 +123,22 @@ qˡ = AquaSkyLES.CondensateField(model, T) qᵛ★ = AquaSkyLES.SaturationField(model, T) δ = Field(model.tracers.q - qᵛ★) +# Output BCs +# Sensible heat flux +Jθ = Oceananigans.Models.BoundaryConditionOperation(model.tracers.θ, :bottom, model) +ρ₀ = AquaSkyLES.base_density(buoyancy) # air density at z=0 +cₚ = buoyancy.thermodynamics.dry_air.heat_capacity +JSH = Field(ρ₀ * cₚ * Jθ) # Heat flux [W/m^2] + +# Latent heat flux +Jq = Oceananigans.Models.BoundaryConditionOperation(model.tracers.q, :bottom, model) +Lᵛ = parameters.condensation.latent_heat +JLH = Field(ρ₀ * Lᵛ * Jq) # [W/m^2] + +# Kinematic momentum flux +Ju = Oceananigans.Models.BoundaryConditionOperation(model.velocities.u, :bottom, model) # [m^2/s^2] +Jv = Oceananigans.Models.BoundaryConditionOperation(model.velocities.v, :bottom, model) # [m^2/s^2] + function progress(sim) compute!(T) @@ -157,7 +173,7 @@ end add_callback!(simulation, progress, IterationInterval(10)) -outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★)) +outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★, JSH, JLH, Ju, Jv)) ow = JLD2Writer(model, outputs, filename = "prescribed_sst_convection.jld2", @@ -177,7 +193,7 @@ Nt = length(θt) using GLMakie, Printf -n = Observable(length(θt)) +n = Observable(1) θn = @lift θt[$n] qn = @lift qt[$n] @@ -200,7 +216,7 @@ hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ)) hmq = heatmap!(axq, qn, colorrange=(0.97e-2, 1.05e-2), colormap=:magma) hmT = heatmap!(axT, Tn, colorrange=(Tmin, Tmax)) hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, 1.5e-3), colormap=:magma) - + # Label(fig[0, 1], "θ", tellwidth=false) # Label(fig[0, 2], "q", tellwidth=false) # Label(fig[0, 1], "θ", tellwidth=false) @@ -216,3 +232,54 @@ fig record(fig, "prescribed_sst.mp4", 1:Nt, framerate=12) do nn n[] = nn end + + +# surface flux timeseries +QSH = FieldTimeSeries("prescribed_sst_convection.jld2", "JSH") +QLH = FieldTimeSeries("prescribed_sst_convection.jld2", "JLH") +Qu = FieldTimeSeries("prescribed_sst_convection.jld2", "Ju") +Qv = FieldTimeSeries("prescribed_sst_convection.jld2", "Jv") + +QSH_av = zeros(length(QSH.times)) +QLH_av = zeros(length(QSH.times)) +ΔT_av = zeros(length(QSH.times)) +Qv_av = zeros(length(QSH.times)) +Qu_av = zeros(length(QSH.times)) + +for n in 1:length(QSH.times) + QSHn = Field(Average(QSH[n], dims=1)) + compute!(QSHn) + QSH_av[n] = QSHn[1, 1, 1] + + QLHn = Field(Average(QLH[n], dims=1)) + compute!(QLHn) + QLH_av[n] = QLHn[1, 1, 1] + + Qun = Field(Average(Qu[n], dims=1)) + compute!(Qun) + Qu_av[n] = Qun[1, 1, 1] + + Qvn = Field(Average(Qv[n], dims=1)) + compute!(Qvn) + Qv_av[n] = Qvn[1, 1, 1] + + ΔT_av[n] = -mean(θt[:,1,1,n])+parameters.sea_surface_temperature + +end + +ΔTtitle = string(znodes(grid, Center())[1], "m air-sea temperature difference") +fig = Figure(size=(600, 700), fontsize=12) +axtau = Axis(fig[1, 1], ylabel=L"\tau/\rho_0 ~(\mathrm{m}^2/\mathrm{s}^2)", xlabel = "time [h]") +axΔT = Axis(fig[3, 1], ylabel=L"\Delta T ~(\mathrm{K})", xlabel = "time [h]", title = ΔTtitle) +axSH = Axis(fig[2, 1], ylabel=L"J~(\mathrm{W}/\mathrm{m}^2)", xlabel = "time [h]") + +lines!(axtau, QSH.times/3600, Qu_av, label = L"\tau_x/\rho_0") +lines!(axtau, QSH.times/3600, Qv_av, label = L"\tau_y/\rho_0") +lines!(axΔT, QSH.times/3600, ΔT_av) +lines!(axSH, QSH.times/3600, QSH_av, label = L"J^{SH}") +lines!(axSH, QLH.times/3600, QLH_av, label = L"J^{LH}") + +fig[1, 2] = Legend(fig, axtau, framevisible = false) +fig[2, 2] = Legend(fig, axSH, framevisible = false) +save("Surface_fluxes.png", fig, px_per_unit = 4) +