Onflyupdates

src/TwoLayerDirectNumericalShenanigans.jl

@@ -1,8 +1,11 @@
 module TwoLayerDirectNumericalShenanigans
 
 using Oceananigans, Printf, Reexport, JLD2, Rasters, NCDatasets, GibbsSeaWater
-using Oceananigans: AbstractModel, Operators.ℑzᵃᵃᶠ
+using Oceananigans: AbstractModel, Operators.ℑzᵃᵃᶜ
+using Oceananigans: BuoyancyModels.get_temperature_and_salinity, BuoyancyModels.θ_and_sᴬ,  BuoyancyModels.Zᶜᶜᶜ
+using SeawaterPolynomials
 using SeawaterPolynomials: TEOS10EquationOfState
+import SeawaterPolynomials.ρ
 using SpecialFunctions: erf
 using Oceanostics: KineticEnergyDissipationRate
 using OceanRasterConversions: get_σₚ

src/kernelfunctions.jl

@@ -1,48 +1,38 @@
+"Extend `ρ′` to compute at user defined reference geopotential height"
+SeawaterPolynomials.ρ(i, j, k, grid, eos, θ, sᴬ) = ρ(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᶜᶜᶜ(i, j, k, grid),eos)
+SeawaterPolynomials.ρ(i, j, k, grid, eos, θ, sᴬ, Zᵣ) = ρ(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᵣ, eos)
 """
-    C(i, j, k, grid, C)
-Get tracer `C` values for use in other function. There may be another way to do this for
-`KernelFunctionOperation`s but I have not found it so will use this for now. **Note:** this
-return the value of the tracer with no interpolation so if the tracer `C` is at
-`(Center, Center, Center)` the value extracted will be at (Center, Center, Center)`.
+    function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
+Compute the density of seawater at grid point `(i, j, k)` using `SeawaterBuoyancy`.
 """
-C(i, j, k, grid, C) = C[i, j, k]
-"""
-    σ(i, j, k, grid, model, reference_pressure)
-Compute potential density `σ` at `reference_pressure` from salinity and temperature tracers
-in `model`.
-"""
-σ(i, j, k, grid, model, reference_pressure) = gsw_rho(C(i, j, k, grid, model.tracers.S),
-                                                      C(i, j, k, grid, model.tracers.T),
-                                                      reference_pressure)
-"""
-    DensityField(model, reference_pressure)
-Return an `KernelFunctionOperation` at `(Center, Center, Center)` that computes the
-potential density from the salinity and temperature tracers in `model` at `reference_pressure`.
-"""
-DensityField(model, reference_pressure) =
-    KernelFunctionOperation{Center, Center, Center}(σ, model.grid, model,
-                                                    reference_pressure)
-"""
-    σ_anomalyᶜᶜᶠ(i, j, k, grid, model, reference_pressure)
-Compute potential density anomaly at `reference_pressure` from salinity and temperature tracers
-in `model` that have been interpolated from `(Center, Center, Center)` to
-`(Center, Center, Face)`.
-"""
-σ_anomalyᶜᶜᶠ(i, j, k, grid, model, reference_pressure, σ_reference) =
-    gsw_rho(ℑzᵃᵃᶠ(i, j, k, grid, model.tracers.S), ℑzᵃᵃᶠ(i, j, k, grid, model.tracers.T),
-            reference_pressure) - σ_reference
-"""
-    function InterpolatedDensityAnomaly(model, reference_pressure)
-Return a `KernelFunctionOperation` at `(Center, Center, Face)` to compute the potential
-density anomaly from the salinity and temperature tracers from `model` that have been
-interpolated from `(Center, Center, Center)` to `(Center, Center, Face)`
-"""
-function InterpolatedDensityAnomaly(model, reference_pressure)
-
-    σ_reference = model.buoyancy.model.equation_of_state.reference_density
-    return KernelFunctionOperation{Center, Center, Face}(σ_anomalyᶜᶜᶠ,
-                                                         model.grid, model,
-                                                         reference_pressure,
-                                                         σ_reference)
-
+@inline function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
+    T, S = get_temperature_and_salinity(b, C)
+    return ρ(i, j, k, grid, b.equation_of_state, T, S)
+end
+density(model) = density(model.buoyancy, model.grid, model.tracers)
+density(b, grid, tracers) =
+    KernelFunctionOperation{Center, Center, Center}(densityᶜᶜᶜ, grid, b.model, tracers)
+Density(model) = density(model)
+"""
+    function potential_densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C, parameters)
+Compute the potential density of seawater at grid point `(i, j, k)`
+at reference pressure `parameters.pᵣ` using `SeawaterBuoyancy`.
+"""
+@inline function potential_densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C, parameters)
+    T, S = get_temperature_and_salinity(b, C)
+    Zᵣ = parameters.Zᵣ
+    return ρ(i, j, k, grid, b.equation_of_state, T, S, Zᵣ)
 end
+potential_density(model, parameters) = potential_density(model.buoyancy, model.grid,
+                                                         model.tracers, parameters)
+potential_density(b, grid, tracers, parameters) =
+    KernelFunctionOperation{Center, Center, Center}(potential_densityᶜᶜᶜ, grid, b.model,
+                                                    tracers, parameters)
+PotentialDensity(model, parameters) = potential_density(model, parameters)
+
+"`(Center, Center, Center)` vertical velocity `Field`"
+wᶜᶜᶜ(model) = wᶜᶜᶜ(model.veolcities.w, model.grid)
+wᶜᶜᶜ(w, grid) = KernelFunctionOperation{Center, Center, Center}(ℑzᵃᵃᶜ, grid, w)
+"`(Center, Center, Face)` vertical buoyancy gradient `Field`"
+∂b∂z(model) = ∂b∂z(model.buoyancy, model.grid, model.tracers)
+∂b∂z(b, grid, tracers) = KernelFunctionOperation{Center, Center, Face}(∂z_b, grid, b, tracers)

src/twolayerdns.jl

@@ -174,31 +174,40 @@ function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number,
 
     # model tracers
     S, T = model.tracers.S, model.tracers.T
-    # custom saved output
 
-    # Density
-    σ = DensityField(model, density_reference_pressure)
+    # Custom saved output
+    # Potential density
+    parameters = (Zᵣ = density_reference_pressure,)
+    σ = PotentialDensity(model, parameters)
 
     # Inferred vertical diffusivity
-    σ_anomaly_interpolated = InterpolatedDensityAnomaly(model, density_reference_pressure)
-    w = model.velocities.w
-    κᵥ = Integral((-w * σ_anomaly_interpolated) / σ)
+    # b_field = BuoyancyField(model)
+    # w_center_field = wᶜᶜᶜ(model)
+    # b_grad_field = ∂b∂z(model)
+    # κᵥ_field = Field((-w_center_field * b_field) / b_grad_field)
+    # compute!(κᵥ_field)
+    # κᵥ = Integral(κᵥ_field)
 
     # Minimum in space Kolmogorov length scale
     ϵ = KineticEnergyDissipationRate(model)
-    η_space(model) = minimum(model.closure.ν ./ ϵ)
+    η_space(model) = (model.closure.ν^3 / maximum(ϵ))^(1/4)
+
+    # Volume integrated TKE dissipation
+    ∫ϵ = Integral(ϵ)
 
     # Dimensions and attributes for custom saved output
-    dims = Dict("η_space" => (), "σ" => ("xC", "xC", "zC"), "κᵥ" => ())
+    dims = Dict("η_space" => (), "σ" => ("xC", "xC", "zC"), #="κᵥ" => (),=# "∫ϵ" => ())
     oa = Dict(
         "σ" => Dict("longname" => "Seawater potential density calculated using TEOS-10 at $(density_reference_pressure)dbar",
                     "units" => "kgm⁻³"),
         "η_space" => Dict("longname" => "Minimum (in space) Kolmogorov length"),
-        "κᵥ" => Dict("longname" => "Inferred vertical diffusivity",
-                     "units" => "m²s⁻¹"))
+        # "κᵥ" => Dict("longname" => "Inferred vertical diffusivity",
+        #              "units" => "m²s⁻¹"),
+        "∫ϵ" => Dict("longname" => "Volume integrated turbulent kintetic energy dissipation")
+        )
 
     # outputs to be saved during the simulation
-    outputs = Dict("S" => S, "T" => T, "η_space" => η_space, "σ" => σ, "κᵥ" => κᵥ)
+    outputs = Dict("S" => S, "T" => T, "η_space" => η_space, "σ" => σ, #="κᵥ" => κᵥ,=# "∫ϵ" => ∫ϵ)
     if save_velocities
         u, v = model.velocities.u, model.velocities.v
         velocities = Dict("u" => u, "v" => v, "w" => w)

test/initialconditions_test.jl

@@ -1,5 +1,4 @@
 # Need to write these to check things are set at the right level in the right field
-architecture = CPU() # or GPU()
 diffusivities = (ν = 1e-4, κ = (S = 1e-5, T = 1e-5))
 resolution = (Nx = 10, Ny = 10, Nz = 500)
 model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;

test/kernelfunction_test.jl

@@ -1,8 +1,9 @@
-using TwoLayerDirectNumericalShenanigans: DensityField, InterpolatedDensityAnomaly
+using TwoLayerDirectNumericalShenanigans: PotentialDensity, Density
 using Oceananigans: Operators.ℑzᵃᵃᶠ
+using Oceananigans: BuoyancyModels.θ_and_sᴬ
+
 using GibbsSeaWater
 
-architecture = CPU()
 diffusivities = (ν = 1e-4, κ = (S = 1e-5, T = 1e-5))
 resolution = (Nx = 10, Ny = 10, Nz = 100)
 
@@ -23,24 +24,26 @@ dns = TwoLayerDNS(model, profile_function, initial_conditions)
 set_two_layer_initial_conditions!(dns)
 
 reference_pressure = 0
-σ_field = Field(DensityField(dns.model, 0))
-compute!(σ_field)
+parameters = (Zᵣ = 0,)
+model = dns.model
+pd_field = Field(PotentialDensity(model, parameters))
+compute!(pd_field)
 
-function test_density_profile(σ_field, computed_density_profile)
+function test_potential_density_profile(pd_field, computed_density_profile, atol)
 
-    vertical_slice = interior(σ_field, rand(1:10), rand(1:10), :)
+    vertical_slice = interior(pd_field, rand(1:10), rand(1:10), :)
 
-    return vertical_slice .== computed_density_profile
+    return isapprox.(vertical_slice, computed_density_profile; atol)
 
 end
 
-σ_anomaly_interp_field = Field(InterpolatedDensityAnomaly(dns.model, reference_pressure))
-compute!(σ_anomaly_interp_field)
+d_field = Field(Density(model))
+compute!(d_field)
 
-function test_density_anom_profile(σ_anomaly_interp_field, computed_density_profile)
+function test_density_profile(d_field, computed_density_profile, atol)
 
-    vertical_slice = interior(σ_anomaly_interp_field, rand(1:10), rand(1:10), :)
+    vertical_slice = interior(d_field, rand(1:10), rand(1:10), :)
 
-    return vertical_slice .== computed_density_profile
+    return isapprox.(vertical_slice, computed_density_profile; atol)
 
 end

test/runtests.jl

@@ -1,6 +1,11 @@
 using TwoLayerDirectNumericalShenanigans, Test
 using TwoLayerDirectNumericalShenanigans: perturb_tracer
+using SeawaterPolynomials
+import SeawaterPolynomials.ρ
+using GibbsSeaWater: gsw_p_from_z
+using CUDA: has_cuda_gpu
 
+architecture = has_cuda_gpu() ? GPU() : CPU()
 include("twolayercontainers_test.jl")
 
 @testset "TwoLayerDNS upper layer containers" begin
@@ -131,34 +136,22 @@ end
 
 include("kernelfunction_test.jl")
 
-@testset "Density field computation" begin
-    d_idx = findfirst(z .== transition_depth)
-    σᵘ = gsw_rho.(interior(model.tracers.S, rand(1:10), rand(1:10), 1:d_idx),
-                  interior(model.tracers.T, rand(1:10), rand(1:10), 1:d_idx),
-                  reference_pressure)
-    σˡ = gsw_rho.(interior(model.tracers.S, rand(1:10), rand(1:10), d_idx+1:100),
-                  interior(model.tracers.T, rand(1:10), rand(1:10), d_idx+1:100),
-                  reference_pressure)
-    σ_computed_profile = vcat(σᵘ, σˡ)
-    @test isequal(trues(length(z)), test_density_profile(σ_field, σ_computed_profile))
+atol = 2e-4 # tolerance for accuracy of density compared to GSW
+@testset "PotentialDensity field computation" begin
+
+    σ_profile = gsw_rho.(interior(model.tracers.S, rand(1:10), rand(1:10), :),
+                         interior(model.tracers.T, rand(1:10), rand(1:10), :),
+                         reference_pressure)
+
+    @test isequal(trues(length(z)),
+                  test_potential_density_profile(pd_field, σ_profile, atol))
 end
 
-@testset "Interpolated density anomaly" begin
-    d_idx = findfirst(z .== transition_depth)
-    S_interpolation = Field(KernelFunctionOperation{Center, Center, Face}(ℑzᵃᵃᶠ, model.grid,
-                                                                          model.tracers.S))
-    compute!(S_interpolation)
-    T_interpolation = Field(KernelFunctionOperation{Center, Center, Face}(ℑzᵃᵃᶠ, model.grid,
-                                                                          model.tracers.T))
-    compute!(T_interpolation)
-    σ_reference = dns.model.buoyancy.model.equation_of_state.reference_density
-    σᵘ = gsw_rho.(interior(S_interpolation, rand(1:10), rand(1:10), 1:d_idx),
-                  interior(T_interpolation, rand(1:10), rand(1:10), 1:d_idx),
-                  reference_pressure) .- σ_reference
-    σˡ = gsw_rho.(interior(S_interpolation, rand(1:10), rand(1:10), d_idx+1:101),
-                  interior(T_interpolation, rand(1:10), rand(1:10), d_idx+1:101),
-                  reference_pressure) .- σ_reference
-    σ_anom_computed_profile = vcat(σᵘ, σˡ)
-    @test isequal(trues(length(σ_anom_computed_profile)),
-                  test_density_anom_profile(σ_anomaly_interp_field, σ_anom_computed_profile))
+@testset "Density field computation" begin
+
+    p = gsw_p_from_z.(z, 0)
+    ρ_profile = gsw_rho.(interior(model.tracers.S, rand(1:10), rand(1:10), :),
+                         interior(model.tracers.T, rand(1:10), rand(1:10), :), p)
+    @test isequal(trues(length(ρ_profile)),
+                   test_density_profile(d_field, ρ_profile, atol))
 end