Merge pull request #140 from jbisits/updateinferredverticaldiffusivity

Change methods for computing IVD

Author: jbisits

Project.toml

@@ -1,7 +1,7 @@
 name = "TwoLayerDirectNumericalShenanigans"
 uuid = "40aaee9f-3595-48be-b36c-f1067009652f"
 authors = ["Josef Bisits <[email protected]>"]
-version = "0.4.3"
+version = "0.4.5"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

src/kernelfunctions.jl

@@ -37,6 +37,17 @@ wᶜᶜᶜ(w, grid) = KernelFunctionOperation{Center, Center, Center}(ℑzᵃᵃ
 ∂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)
 
+@inline vertical_buoyancy_flux(i, j, k, grid, b::SeawaterBuoyancy, C, w) =
+        -ℑzᵃᵃᶜ(i, j, k, w) * buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, b, C)
+@inline function vertical_buoyancy_flux(model)
+
+    grid = model.grid
+    b = model.buoyancy.model
+    C = model.tracers
+    w = model.velocities.w
+
+    return KernelFunctionOperation{Center, Center, Center}(vertical_buoyancy_flux, grid, b, C, w)
+end
 @inline function Kᵥ(i, j, k, grid, b::SeawaterBuoyancy, C, w)
 
     if ∂z_b(i, j, k, grid, b, C) != 0

src/output_utils.jl

@@ -331,6 +331,39 @@ function non_dimensional_numbers!(simulation::Simulation, dns::TwoLayerDNS)
 
     return nothing
 
+end
+"""
+    function inferred_vertical_diffusivity(saved_output::AbstractString)
+Calculate the inferred vertical diffusivity using the horizontally averaged vertical
+buoyancy gradient and horizontally averaged vertical buoyancy flux
+```math
+κᵥ = -\\frac{\\overline{w'b'}}{\\overline{\\frac{∂b}{∂z}}}.
+```
+**Note:** the horizontally avergaed vertical buoyancy gradient and 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)
+
+    NCDataset(saved_output, "a") do ds
+        b_grad = ds[:∫ₐb_grad][2:end, :]
+        replace!(b_grad, 0 => NaN)
+        b_flux = ds[:∫ₐb_flux][:, :]
+        ∫ₐκᵥ = similar(b_flux)
+        ∫ₐκᵥ .= b_flux ./ b_grad
+        defVar(ds, "∫ₐκᵥ", ∫ₐκᵥ, ("zC", "time"),
+               attrib = Dict("longname" => "Horizontally integrated inferred vertical diffusivity",
+                             "units" => "m²s⁻¹"))
+        dV = (diff(ds[:xC][1:2]) .* diff(ds[:yC][1:2])) .* diff(ds[:zF][:])
+        replace!(∫ₐκᵥ, NaN => 0)
+        ∫κᵥ = mapslices(sum, ∫ₐκᵥ .* dV, dims = 1)
+        defVar(ds, "∫κᵥ", ∫κᵥ, ("time",),
+        attrib = Dict("longname" => "Volume integrated inferred vertical diffusivity",
+                      "units" => "m²s⁻¹"))
+    end
+
+    return nothing
+
 end
 """
     funciton find_file_type(file::AbstractString)

src/twolayerdns.jl

@@ -181,29 +181,34 @@ function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number,
     σ = PotentialDensity(model, parameters)
 
     # Inferred vertical diffusivity
-    κᵥ = InferredVerticalDiffusivity(model)
-    ∫κᵥ = Integral(κᵥ)
+    b_flux = vertical_buoyancy_flux(model)
+    ∫ₐb_flux = Integral(b_flux, dims = (1, 2))
+    b_grad = ∂b∂z(model)
+    ∫ₐb_grad = Integral(b_grad, dims = (1, 2))
 
-    # Minimum in space Kolmogorov length scale
     ϵ = KineticEnergyDissipationRate(model)
+    # Volume integrated TKE dissipation
+    ∫ϵ = Integral(ϵ)
+    # Minimum in space Kolmogorov length scale
     η_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" => ())
     oa = Dict(
         "σ" => Dict("longname" => "Seawater potential density calculated using TEOS-10 at $(density_reference_gp_height)dbar",
                     "units" => "kgm⁻³"),
         "η_space" => Dict("longname" => "Minimum (in space) Kolmogorov length"),
-        "∫κᵥ" => Dict("longname" => "Volume integrated inferred vertical diffusivity",
-                      "units" => "m²s⁻¹"),
+        "∫ₐb_flux" => Dict("longname" => "Horizontally integrated vertical buoyacny flux (w'b')"),
+        "∫ₐb_grad" => Dict("longname" => "Horizontally integrated vertical buoyancy gradient (∂b/∂z)"),
         "∫ϵ" => 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, "σ" => σ, "∫ϵ" => ∫ϵ,
+                  "∫ₐb_flux" => ∫ₐb_flux, "∫ₐb_grad" => ∫ₐb_grad)
     if save_velocities
         u, v, w = model.velocities
         velocities = Dict("u" => u, "v" => v, "w" => w)