From e54eb05e1fe77ca0ea1dbc558c225f61d7744e98 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Sun, 23 Feb 2025 14:47:50 -0800 Subject: [PATCH 1/4] first test with PV --- src/FlowDiagnostics.jl | 26 +++++++++++++++++++++++++- src/Oceanostics.jl | 4 ++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/src/FlowDiagnostics.jl b/src/FlowDiagnostics.jl index 0b2ee3a1..aacdcbd5 100755 --- a/src/FlowDiagnostics.jl +++ b/src/FlowDiagnostics.jl @@ -5,7 +5,8 @@ export RichardsonNumber, RossbyNumber export ErtelPotentialVorticity, ThermalWindPotentialVorticity, DirectionalErtelPotentialVorticity export StrainRateTensorModulus, VorticityTensorModulus, Q, QVelocityGradientTensorInvariant -using Oceanostics: validate_location, +using Oceanostics: AbstractDiagnostic, + validate_location, validate_dissipative_closure, add_background_fields, get_coriolis_frequency_components @@ -283,6 +284,29 @@ function ErtelPotentialVorticity(model, u, v, w, tracer, coriolis; location = (F u, v, w, tracer, fx, fy, fz) end +struct PotentialVorticity <: AbstractDiagnostic + operation + thermal_wind + tracer +end + +Base.summary(pv::PotentialVorticity) = string("PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with $(summary(pv.operation))") +function Base.show(io::IO, pv::PotentialVorticity) + print(io, "PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with\n") + show(io, pv.operation) +end + + +function PotentialVorticity(model, u, v, w, tracer, coriolis; thermal_wind = false, location = (Face, Face, Face)) + validate_location(location, "PotentialVorticity", (Face, Face, Face)) + if thermal_wind + kfo = ThermalWindPotentialVorticity(model, u, v, tracer, coriolis; location) + else + kfo = ErtelPotentialVorticity(model, u, v, w, tracer, coriolis; location) + end + return PotentialVorticity(kfo, thermal_wind, tracer) +end + @inline function directional_ertel_potential_vorticity_fff(i, j, k, grid, u, v, w, b, params) dWdy = ℑxᶠᵃᵃ(i, j, k, grid, ∂yᶜᶠᶠ, w) # C, C, F → C, F, F → F, F, F diff --git a/src/Oceanostics.jl b/src/Oceanostics.jl index 86238c0c..453c6562 100755 --- a/src/Oceanostics.jl +++ b/src/Oceanostics.jl @@ -30,6 +30,10 @@ export PotentialEnergy export ProgressMessengers #--- +#+++ Define AbstractDiagnostic +abstract type AbstractDiagnostic end +#--- + #+++ Utils for validation # Right now, all kernels must be located at ccc using Oceananigans.TurbulenceClosures: AbstractScalarDiffusivity, ThreeDimensionalFormulation From aef0d4ef8447f8e7b40f076512893c06a91c6632 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Mon, 24 Feb 2025 19:32:43 -0800 Subject: [PATCH 2/4] more methods --- src/FlowDiagnostics.jl | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/FlowDiagnostics.jl b/src/FlowDiagnostics.jl index aacdcbd5..3c91b6d5 100755 --- a/src/FlowDiagnostics.jl +++ b/src/FlowDiagnostics.jl @@ -290,23 +290,29 @@ struct PotentialVorticity <: AbstractDiagnostic tracer end -Base.summary(pv::PotentialVorticity) = string("PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with $(summary(pv.operation))") -function Base.show(io::IO, pv::PotentialVorticity) - print(io, "PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with\n") - show(io, pv.operation) -end - +PotentialVorticity(model; tracer, location = (Face, Face, Face)) = PotentialVorticity(model, model.velocities..., tracer, model.coriolis; location) +PotentialVorticity(model, tracer::Symbol) = PotentialVorticity(model, model.velocities..., model.tracers[tracer], model.coriolis; location) function PotentialVorticity(model, u, v, w, tracer, coriolis; thermal_wind = false, location = (Face, Face, Face)) validate_location(location, "PotentialVorticity", (Face, Face, Face)) + f⃗ = get_coriolis_frequency_components(coriolis) if thermal_wind - kfo = ThermalWindPotentialVorticity(model, u, v, tracer, coriolis; location) + return KernelFunctionOperation{Face, Face, Face}(potential_vorticity_in_thermal_wind_fff, model.grid, + u, v, tracer, f⃗...) else - kfo = ErtelPotentialVorticity(model, u, v, w, tracer, coriolis; location) + return KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_fff, model.grid, + u, v, w, tracer, f⃗...) end return PotentialVorticity(kfo, thermal_wind, tracer) end +Base.summary(pv::PotentialVorticity) = string("PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with $(summary(pv.operation))") +function Base.show(io::IO, pv::PotentialVorticity) + print(io, "PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with\n") + show(io, pv.operation) +end + + @inline function directional_ertel_potential_vorticity_fff(i, j, k, grid, u, v, w, b, params) dWdy = ℑxᶠᵃᵃ(i, j, k, grid, ∂yᶜᶠᶠ, w) # C, C, F → C, F, F → F, F, F From 988ab76b122b5aeff4e8cfea35c450ab3eba92b6 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Mon, 24 Feb 2025 19:46:16 -0800 Subject: [PATCH 3/4] go back to Ertel potential vorticity --- src/FlowDiagnostics.jl | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/src/FlowDiagnostics.jl b/src/FlowDiagnostics.jl index 3c91b6d5..d558a5e1 100755 --- a/src/FlowDiagnostics.jl +++ b/src/FlowDiagnostics.jl @@ -272,29 +272,17 @@ Note that EPV values are correctly calculated both in the interior and the bound interior and top boundary, EPV = f×N² = 10⁻¹⁰, while EPV = 0 at the bottom boundary since ∂b/∂z is zero there. """ -function ErtelPotentialVorticity(model; tracer = :b, location = (Face, Face, Face)) - validate_location(location, "ErtelPotentialVorticity", (Face, Face, Face)) - return ErtelPotentialVorticity(model, model.velocities..., model.tracers[tracer], model.coriolis; location) -end - -function ErtelPotentialVorticity(model, u, v, w, tracer, coriolis; location = (Face, Face, Face)) - validate_location(location, "ErtelPotentialVorticity", (Face, Face, Face)) - fx, fy, fz = get_coriolis_frequency_components(coriolis) - return KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_fff, model.grid, - u, v, w, tracer, fx, fy, fz) -end - -struct PotentialVorticity <: AbstractDiagnostic +struct ErtelPotentialVorticity <: AbstractDiagnostic operation thermal_wind tracer end -PotentialVorticity(model; tracer, location = (Face, Face, Face)) = PotentialVorticity(model, model.velocities..., tracer, model.coriolis; location) -PotentialVorticity(model, tracer::Symbol) = PotentialVorticity(model, model.velocities..., model.tracers[tracer], model.coriolis; location) +ErtelPotentialVorticity(model; tracer::Symbol, location = (Face, Face, Face)) = ErtelPotentialVorticity(model, model.velocities..., model.tracers[tracer], model.coriolis; location) +ErtelPotentialVorticity(model, tracer; location = (Face, Face, Face)) = ErtelPotentialVorticity(model, model.velocities..., tracer, model.coriolis; location) -function PotentialVorticity(model, u, v, w, tracer, coriolis; thermal_wind = false, location = (Face, Face, Face)) - validate_location(location, "PotentialVorticity", (Face, Face, Face)) +function ErtelPotentialVorticity(model, u, v, w, tracer, coriolis; thermal_wind = false, location = (Face, Face, Face)) + validate_location(location, "ErtelPotentialVorticity", (Face, Face, Face)) f⃗ = get_coriolis_frequency_components(coriolis) if thermal_wind return KernelFunctionOperation{Face, Face, Face}(potential_vorticity_in_thermal_wind_fff, model.grid, @@ -303,12 +291,12 @@ function PotentialVorticity(model, u, v, w, tracer, coriolis; thermal_wind = fal return KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_fff, model.grid, u, v, w, tracer, f⃗...) end - return PotentialVorticity(kfo, thermal_wind, tracer) + return ErtelPotentialVorticity(kfo, thermal_wind, tracer) end -Base.summary(pv::PotentialVorticity) = string("PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with $(summary(pv.operation))") -function Base.show(io::IO, pv::PotentialVorticity) - print(io, "PotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with\n") +Base.summary(pv::ErtelPotentialVorticity) = string("ErtelPotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with $(summary(pv.operation))") +function Base.show(io::IO, pv::ErtelPotentialVorticity) + print(io, "ErtelPotentialVorticity$(pv.thermal_wind ? " in thermal wind balance" : ""), calculated with\n") show(io, pv.operation) end From 0d184349fba73d441defef7d86ea2142e22bba41 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Tue, 25 Feb 2025 20:09:16 -0800 Subject: [PATCH 4/4] fix bug and merge methods --- src/FlowDiagnostics.jl | 48 +++++++++--------------------------------- 1 file changed, 10 insertions(+), 38 deletions(-) diff --git a/src/FlowDiagnostics.jl b/src/FlowDiagnostics.jl index d558a5e1..7b1495ab 100755 --- a/src/FlowDiagnostics.jl +++ b/src/FlowDiagnostics.jl @@ -154,8 +154,7 @@ end #--- #+++ Potential vorticity -@inline function potential_vorticity_in_thermal_wind_fff(i, j, k, grid, u, v, b, f) - +@inline function ertel_potential_vorticity_in_thermal_wind_fff(i, j, k, grid, u, v, b, f) dVdx = ℑzᵃᵃᶠ(i, j, k, grid, ∂xᶠᶠᶜ, v) # F, F, C → F, F, F dUdy = ℑzᵃᵃᶠ(i, j, k, grid, ∂yᶠᶠᶜ, u) # F, F, C → F, F, F dbdz = ℑxyᶠᶠᵃ(i, j, k, grid, ∂zᶜᶜᶠ, b) # C, C, F → F, F, F @@ -170,32 +169,6 @@ end return pv_barot + pv_baroc end -""" - $(SIGNATURES) - -Calculate the Potential Vorticty assuming thermal wind balance for `model`, where the characteristics of -the Coriolis rotation are taken from `model.coriolis`. The Potential Vorticity in this case -is defined as - -``` - TWPV = (f + ωᶻ) ∂b/∂z - f ((∂U/∂z)² + (∂V/∂z)²) -``` -where `f` is the Coriolis frequency, `ωᶻ` is the relative vorticity in the `z` direction, `b` is the buoyancy, and -`∂U/∂z` and `∂V/∂z` comprise the thermal wind shear. -""" -function ThermalWindPotentialVorticity(model; tracer = :b, location = (Face, Face, Face)) - validate_location(location, "ThermalWindPotentialVorticity", (Face, Face, Face)) - u, v, w = model.velocities - return ThermalWindPotentialVorticity(model, u, v, model.tracers[tracer], model.coriolis; location) -end - -function ThermalWindPotentialVorticity(model, u, v, tracer, coriolis; location = (Face, Face, Face)) - validate_location(location, "ThermalWindPotentialVorticity", (Face, Face, Face)) - fx, fy, fz = get_coriolis_frequency_components(coriolis) - return KernelFunctionOperation{Face, Face, Face}(potential_vorticity_in_thermal_wind_fff, model.grid, - u, v, tracer, fz) -end - @inline function ertel_potential_vorticity_fff(i, j, k, grid, u, v, w, b, fx, fy, fz) dWdy = ℑxᶠᵃᵃ(i, j, k, grid, ∂yᶜᶠᶠ, w) # C, C, F → C, F, F → F, F, F dVdz = ℑxᶠᵃᵃ(i, j, k, grid, ∂zᶜᶠᶠ, v) # C, F, C → C, F, F → F, F, F @@ -218,14 +191,13 @@ end """ $(SIGNATURES) -Calculate the Ertel Potential Vorticty for `model`, where the characteristics of -the Coriolis rotation are taken from `model.coriolis`. The Ertel Potential Vorticity -is defined as - EPV = ωₜₒₜ ⋅ ∇b +Calculate the Ertel Potential Vorticty, defined as -where ωₜₒₜ is the total (relative + planetary) vorticity vector, `b` is the buoyancy and ∇ is the gradient -operator. + EPV = ωₜₒₜ ⋅ ∇c + +(where ωₜₒₜ is the total (relative + planetary) vorticity vector, `c` is a conserved tracer and ∇ is the gradient +operator), for a `model`. ```jldoctest julia> using Oceananigans @@ -285,11 +257,11 @@ function ErtelPotentialVorticity(model, u, v, w, tracer, coriolis; thermal_wind validate_location(location, "ErtelPotentialVorticity", (Face, Face, Face)) f⃗ = get_coriolis_frequency_components(coriolis) if thermal_wind - return KernelFunctionOperation{Face, Face, Face}(potential_vorticity_in_thermal_wind_fff, model.grid, - u, v, tracer, f⃗...) + kfo = KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_in_thermal_wind_fff, model.grid, + u, v, tracer, f⃗...) else - return KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_fff, model.grid, - u, v, w, tracer, f⃗...) + kfo = KernelFunctionOperation{Face, Face, Face}(ertel_potential_vorticity_fff, model.grid, + u, v, w, tracer, f⃗...) end return ErtelPotentialVorticity(kfo, thermal_wind, tracer) end