Missing dep + add ! + use seawater_density

Project.toml

@@ -31,7 +31,7 @@ GibbsSeaWater = "0.1"
 JLD2 = "0.4"
 NCDatasets = "0.12"
 OceanRasterConversions = "0.5"
-Oceananigans = "0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89"
+Oceananigans = "0.89"
 Oceanostics = "0.11, 0.12, 0.13"
 RasterHistograms = "0.1"
 Rasters = "0.8"

ext/TLDNSMakieExt/TLDNSMakieExt.jl

@@ -1,6 +1,6 @@
 module TLDNSMakieExt
 
-using CairoMakie, TwoLayerDirectNumericalShenanigans, Printf,
+using CairoMakie, TwoLayerDirectNumericalShenanigans, Printf, NCDatasets,
       GibbsSeaWater, Rasters, RasterHistograms
 using Oceananigans.Fields
 

ext/TLDNSMakieExt/plotrecipes.jl

@@ -148,6 +148,40 @@ function TLDNS.plot_scalar_diagnostics(saved_output::AbstractString)
 
     return nothing
 
+end
+"""
+    function hovmoller(saved_output::AbstractString, variable::AbstractString)
+Produce a Hovmoller plot of a `variable` in `saved_ouput`.
+"""
+function TLDNS.hovmoller(saved_output::AbstractString, variable::AbstractString;
+                         colormap = :viridis, unit = nothing)
+
+    NCDataset(saved_output) do ds
+
+        t = ds["time"][:] ./ 60
+        zC = ds["zC"][:]
+        plot_var = ds[variable][:, :]
+        fig = Figure(size = (600, 1000))
+        ax = Axis(fig[1, 1],
+                  title = "Hovmoller plot of $(variable)",
+                  xlabel = "t (minutes)",
+                  ylabel = "z (m)")
+        colormap = cgrad(colormap)[2:end-1]
+        lowclip = cgrad(colormap)[1]
+        highclip = cgrad(colormap)[end]
+        hm = heatmap!(ax, t, zC, plot_var';
+                      colormap, lowclip, highclip, colorrange = (-1e-4, 1e-4))
+        cbar_label = isnothing(unit) ? variable : variable * unit
+        Colorbar(fig[1, 2], hm, label = cbar_label)
+
+        plotsave = "hovmoller_"*variable*".png"
+        save(plotsave, fig)
+        @info "Save plot to $(plotsave)"
+
+    end
+
+    return nothing
+
 end
 """
     function animate_volume_distributions(rs::Raster)

src/TwoLayerDirectNumericalShenanigans.jl

@@ -4,6 +4,7 @@ using Oceananigans, Printf, Reexport, JLD2, Rasters, NCDatasets, GibbsSeaWater
 using Oceananigans: AbstractModel, Operators.ℑzᵃᵃᶜ
 using Oceananigans: BuoyancyModels.get_temperature_and_salinity, BuoyancyModels.θ_and_sᴬ,  BuoyancyModels.Zᶜᶜᶜ
 using Oceananigans: BuoyancyModels.buoyancy_perturbationᶜᶜᶜ, BuoyancyModels.∂z_b
+using Oceananigans: Models.seawater_density
 using SeawaterPolynomials
 using SeawaterPolynomials: TEOS10EquationOfState
 import SeawaterPolynomials.ρ
@@ -59,7 +60,7 @@ export DOMAIN_EXTENT, HIGH_RESOLUTION, SO_DIFFUSIVITIES, REFERENCE_DENSITY,
 export compute_density, compute_density!
 
 export animate_2D_field, visualise_initial_conditions, visualise_initial_density,
-       visualise_snapshot, plot_scalar_diagnostics
+       visualise_snapshot, plot_scalar_diagnostics, hovmoller
 
 include("initialconditions.jl")
 include("continuousprofilefunctions.jl")

src/kernelfunctions.jl

@@ -1,35 +1,3 @@
-"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)
-"""
-    function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
-Compute the density of seawater at grid point `(i, j, k)` using `SeawaterBuoyancy`.
-"""
-@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)
@@ -66,3 +34,38 @@ function InferredVerticalDiffusivity(model)
 
     return KernelFunctionOperation{Center, Center, Center}(Kᵥ, grid, b, C, w)
 end
+
+## This has been implemented (by me) in Oceananigans.jl as of v0.89.3. Once I know
+# everything is working I will remove this in favour of Oceananigans.jl version.
+
+# "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)
+# """
+#     function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
+# Compute the density of seawater at grid point `(i, j, k)` using `SeawaterBuoyancy`.
+# """
+# @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)

src/makiefunctions.jl

@@ -8,3 +8,4 @@ function visualise_snapshot end
 function animate_volume_distributions end
 function volume_distribution_snaphsots end
 function plot_scalar_diagnostics end
+function hovmoller end

src/output_utils.jl

@@ -343,7 +343,7 @@ buoyancy gradient and horizontally averaged vertical buoyancy flux
 must be saved in `saved_output`. This is the default behaviour as of version 0.4.5.
 Further it is assumed the horizontal resolution is equal.
 """
-function inferred_vertical_diffusivity(saved_output::AbstractString)
+function inferred_vertical_diffusivity!(saved_output::AbstractString)
 
     NCDataset(saved_output, "a") do ds
         b_grad = ds[:∫ₐb_grad][2:end, :]

src/twolayerdns.jl

@@ -177,8 +177,7 @@ function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number,
 
     # Custom saved output
     # Potential density
-    parameters = (Zᵣ = density_reference_gp_height,)
-    σ = PotentialDensity(model, parameters)
+    σ = seawater_density(model, geopotential_height = 0)
 
     # Inferred vertical diffusivity
     b_flux = vertical_buoyancy_flux(model)

test/runtests.jl

@@ -134,24 +134,24 @@ end
     @test isequal(z[10:20], find_depth(model, test_depth_range))
 end
 
-include("kernelfunction_test.jl")
+# include("kernelfunction_test.jl")
 
-atol = 2e-4 # tolerance for accuracy of density compared to GSW
-@testset "PotentialDensity field computation" begin
+# 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)
+#     σ_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
+#     @test isequal(trues(length(z)),
+#                   test_potential_density_profile(pd_field, σ_profile, atol))
+# end
 
-@testset "Density field computation" begin
+# @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
+#     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