Skip to content

Commit 9df02ef

Browse files
jbisitsjbisits
authored
Merge pull request #160 from jbisits/tracerinferredvdiff
Add `vertical_tracer_flux` + tests
2 parents c02d159 + df60fd8 commit 9df02ef

File tree

6 files changed

+37
-138
lines changed

6 files changed

+37
-138
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "TwoLayerDirectNumericalShenanigans"
22
uuid = "40aaee9f-3595-48be-b36c-f1067009652f"
33
authors = ["Josef Bisits <[email protected]>"]
4-
version = "0.5.1"
4+
version = "0.5.2"
55

66
[deps]
77
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

‎examples/twolayer_example.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@ stop_time = 60
3131
save_schedule = 0.5 # seconds
3232
output_path = joinpath(@__DIR__, "outputs/")
3333
checkpointer_time_interval = 30
34-
simulation = TLDNS_simulation_setup(tldns, Δt, stop_time, save_schedule;
34+
simulation = TLDNS_simulation_setup(tldns, Δt, stop_time, save_schedule, TLDNS.save_computed_output!;
3535
output_path, checkpointer_time_interval)
3636

3737
## Run the simulation

‎src/kernelfunctions.jl‎

Lines changed: 9 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -1,83 +1,21 @@
1-
"`(Center, Center, Center)` vertical velocity `Field`"
2-
wᶜᶜᶜ(model) = wᶜᶜᶜ(model.velocities.w, model.grid)
3-
wᶜᶜᶜ(w, grid) = KernelFunctionOperation{Center, Center, Center}(ℑzᵃᵃᶜ, grid, w)
4-
"`(Center, Center, Face)` vertical buoyancy gradient `Field`"
5-
∂b∂z(model) = ∂b∂z(model.buoyancy, model.grid, model.tracers)
6-
∂b∂z(b, grid, tracers) = KernelFunctionOperation{Center, Center, Face}(∂z_b, grid, b, tracers)
7-
8-
@inline vertical_buoyancy_flux(i, j, k, grid, b::SeawaterBuoyancy, C, w) =
9-
-ℑzᵃᵃᶜ(i, j, k, w) * buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, b, C)
10-
@inline function vertical_buoyancy_flux(model)
11-
12-
grid = model.grid
13-
b = model.buoyancy.model
14-
C = model.tracers
15-
w = model.velocities.w
16-
17-
return KernelFunctionOperation{Center, Center, Center}(vertical_buoyancy_flux, grid, b, C, w)
18-
end
19-
@inline function Kᵥ(i, j, k, grid, b::SeawaterBuoyancy, C, w)
20-
21-
if ∂z_b(i, j, k, grid, b, C) != 0
22-
(-ℑzᵃᵃᶜ(i, j, k, w) * buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, b, C)) / ∂z_b(i, j, k, grid, b, C)
23-
else
24-
0
25-
end
26-
27-
end
28-
function InferredVerticalDiffusivity(model)
1+
@inline tracer_perturbationᶜᶜᶜ(i, j, k, grid, tracer, tracer_mean) =
2+
tracer[i, j, k] - tracer_mean[i, j, k]
3+
# Only needed for testing
4+
@inline function tracer_perturbation(model, tracer)
295

306
grid = model.grid
31-
b = model.buoyancy.model
32-
C = model.tracers
33-
w = model.velocities.w
7+
tracer_mean = Field(Average(tracer))
348

35-
return KernelFunctionOperation{Center, Center, Center}(Kᵥ, grid, b, C, w)
9+
return KernelFunctionOperation{Center, Center, Center}(tracer_perturbationᶜᶜᶜ, grid, tracer, tracer_mean)
3610
end
3711

38-
@inline tracer_perturbationᶜᶜᶜ(i, j, k, grid, tracer, tracer_mean) = tracer[i, j, k] - tracer_mean
39-
@inline vertical_tracer_flux(i, j, k, grid, tracer, tracer_mean, w) =
12+
@inline vtf(i, j, k, grid, tracer, tracer_mean, w) =
4013
-ℑzᵃᵃᶜ(i, j, k, w) * tracer_perturbationᶜᶜᶜ(i, j, k, grid, tracer, tracer_mean)
4114
@inline function vertical_tracer_flux(model, tracer)
4215

4316
grid = model.grid
44-
tracer_mean = Average(tracer)
17+
tracer_mean = Field(Average(tracer))
4518
w = model.velocities.w
4619

47-
return KernelFunctionOperation{Center, Center, Center}(vertical_tracer_flux, grid, tracer, tracer_mean, w)
20+
return KernelFunctionOperation{Center, Center, Center}(vtf, grid, tracer, tracer_mean, w)
4821
end
49-
50-
## This has been implemented (by me) in Oceananigans.jl as of v0.89.3. Once I know
51-
# everything is working I will remove this in favour of Oceananigans.jl version.
52-
53-
# "Extend `ρ′` to compute at user defined reference geopotential height"
54-
# SeawaterPolynomials.ρ(i, j, k, grid, eos, θ, sᴬ) = ρ(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᶜᶜᶜ(i, j, k, grid),eos)
55-
# SeawaterPolynomials.ρ(i, j, k, grid, eos, θ, sᴬ, Zᵣ) = ρ(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᵣ, eos)
56-
# """
57-
# function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
58-
# Compute the density of seawater at grid point `(i, j, k)` using `SeawaterBuoyancy`.
59-
# """
60-
# @inline function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
61-
# T, S = get_temperature_and_salinity(b, C)
62-
# return ρ(i, j, k, grid, b.equation_of_state, T, S)
63-
# end
64-
# density(model) = density(model.buoyancy, model.grid, model.tracers)
65-
# density(b, grid, tracers) =
66-
# KernelFunctionOperation{Center, Center, Center}(densityᶜᶜᶜ, grid, b.model, tracers)
67-
# Density(model) = density(model)
68-
# """
69-
# function potential_densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C, parameters)
70-
# Compute the potential density of seawater at grid point `(i, j, k)`
71-
# at reference pressure `parameters.pᵣ` using `SeawaterBuoyancy`.
72-
# """
73-
# @inline function potential_densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C, parameters)
74-
# T, S = get_temperature_and_salinity(b, C)
75-
# Zᵣ = parameters.Zᵣ
76-
# return ρ(i, j, k, grid, b.equation_of_state, T, S, Zᵣ)
77-
# end
78-
# potential_density(model, parameters) = potential_density(model.buoyancy, model.grid,
79-
# model.tracers, parameters)
80-
# potential_density(b, grid, tracers, parameters) =
81-
# KernelFunctionOperation{Center, Center, Center}(potential_densityᶜᶜᶜ, grid, b.model,
82-
# tracers, parameters)
83-
# PotentialDensity(model, parameters) = potential_density(model, parameters)

‎src/twolayerdns.jl‎

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -307,23 +307,24 @@ function save_computed_output!(simulation, tldns, save_schedule, save_file, outp
307307
ϵ = KineticEnergyDissipationRate(model)
308308
∫ϵ = Integral(ϵ)
309309
ϵ_maximum = Reduction(maximum!, ϵ, dims = (1, 2, 3))
310-
computed_outputs = Dict("σ" => σ, "∫ϵ" => ∫ϵ, "ϵ_maximum" => ϵ_maximum)
310+
311+
# Horizontally integrated vertical temperature flux and vertical temperature gradient
312+
T = model.tracers.T
313+
w′T′ = vertical_tracer_flux(model, T)
314+
∫ₐw′T′ = Integral(w′T′, dims = (1, 2))
315+
∂T∂z = ∂z(T)
316+
∫ₐ∂T∂z = Integral(∂T∂z, dims = (1, 2))
317+
318+
computed_outputs = Dict("σ" => σ, "∫ϵ" => ∫ϵ, "ϵ_maximum" => ϵ_maximum, "∫ₐ∂T∂z" => ∫ₐ∂T∂z,
319+
"∫ₐw′T′" => ∫ₐw′T′)
311320
oa = Dict(
312321
"σ" => Dict("longname" => "Seawater potential density calculated using TEOS-10 at $(reference_gp_height)dbar",
313322
"units" => "kgm⁻³"),
314323
"ϵ_maximum" => Dict("longname" => "Maximum (in space) TKE dissipation"),
315-
"∫ϵ" => Dict("longname" => "Volume integrated turbulent kintetic energy dissipation")
324+
"∫ϵ" => Dict("longname" => "Volume integrated turbulent kintetic energy dissipation"),
325+
"∫ₐw′T′" => Dict("longname" => "Horizontally integrated vertical temeprature flux"),
326+
"∫ₐ∂T∂z" => Dict("longname" => "Horizontally integrated vertical temperature gradient")
316327
)
317-
318-
#TODO: finish this calculation.
319-
# Inferred vertical temperature diffusivity
320-
# T_mean = Average(T)
321-
# T′ = T - Field(T_mean)
322-
# w′ = wᶜᶜᶜ(model)
323-
# w′T′ = Field(w′ * T′)
324-
# ∫ₐw′T′ = Integral(w′T′, dims = (1, 2))
325-
# T_gradient = ∂z(T)
326-
# ∫ₐT_gradient = Integral(T_gradient, dims = (1, 2))
327328
simulation.output_writers[:computed_output] =
328329
save_file == :netcdf ? NetCDFOutputWriter(model, computed_outputs;
329330
filename = "computed_output",

‎test/kernelfunction_test.jl‎

Lines changed: 2 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -1,49 +1,16 @@
1-
using TwoLayerDirectNumericalShenanigans: PotentialDensity, Density
2-
using Oceananigans: Operators.ℑzᵃᵃᶠ
3-
using Oceananigans: BuoyancyModels.θ_and_sᴬ
4-
5-
using GibbsSeaWater
6-
71
diffusivities = (ν = 1e-4, κ = (S = 1e-5, T = 1e-5))
82
resolution = (Nx = 10, Ny = 10, Nz = 100)
93

104
model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
115
reference_density = REFERENCE_DENSITY)
126

13-
z = znodes(model.grid, Center())
14-
157
T₀ᵘ = -1.5
168
S₀ᵘ = 34.568
179
cabbeling = StableUpperLayerInitialConditions(S₀ᵘ, T₀ᵘ)
1810
initial_conditions = TwoLayerInitialConditions(cabbeling)
1911
transition_depth = find_depth(model, INTERFACE_LOCATION)
2012
profile_function = StepChange(transition_depth)
2113

22-
dns = TwoLayerDNS(model, profile_function, initial_conditions)
23-
24-
set_two_layer_initial_conditions!(dns)
25-
26-
reference_pressure = 0
27-
parameters = (Zᵣ = 0,)
28-
model = dns.model
29-
pd_field = Field(PotentialDensity(model, parameters))
30-
compute!(pd_field)
31-
32-
function test_potential_density_profile(pd_field, computed_density_profile, atol)
33-
34-
vertical_slice = interior(pd_field, rand(1:10), rand(1:10), :)
35-
36-
return isapprox.(vertical_slice, computed_density_profile; atol)
37-
38-
end
39-
40-
d_field = Field(Density(model))
41-
compute!(d_field)
42-
43-
function test_density_profile(d_field, computed_density_profile, atol)
44-
45-
vertical_slice = interior(d_field, rand(1:10), rand(1:10), :)
46-
47-
return isapprox.(vertical_slice, computed_density_profile; atol)
14+
tldns = TwoLayerDNS(model, profile_function, initial_conditions)
4815

49-
end
16+
set_two_layer_initial_conditions!(tldns)

‎test/runtests.jl‎

Lines changed: 11 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -180,24 +180,17 @@ include("output_test.jl")
180180

181181
end
182182

183-
# include("kernelfunction_test.jl")
183+
include("kernelfunction_test.jl")
184184

185-
# atol = 2e-4 # tolerance for accuracy of density compared to GSW
186-
# @testset "PotentialDensity field computation" begin
185+
@testset "Vertical temperature flux" begin
187186

188-
# σ_profile = gsw_rho.(interior(model.tracers.S, rand(1:10), rand(1:10), :),
189-
# interior(model.tracers.T, rand(1:10), rand(1:10), :),
190-
# reference_pressure)
187+
T_anom = Field(TLDNS.tracer_perturbation(tldns.model, tldns.model.tracers.T))
188+
compute!(T_anom)
189+
computed_mean = sum(interior(model.tracers.T)) / *(size(interior(model.tracers.T))...)
190+
computed_T_anom = interior(model.tracers.T) .- computed_mean
191+
@test all(computed_T_anom .== interior(T_anom))
192+
vtflux = Field(TLDNS.vertical_tracer_flux(tldns.model, tldns.model.tracers.T))
193+
compute!(vtflux)
194+
@test all(vtflux.data .≈ 0 )
191195

192-
# @test isequal(trues(length(z)),
193-
# test_potential_density_profile(pd_field, σ_profile, atol))
194-
# end
195-
196-
# @testset "Density field computation" begin
197-
198-
# p = gsw_p_from_z.(z, 0)
199-
# ρ_profile = gsw_rho.(interior(model.tracers.S, rand(1:10), rand(1:10), :),
200-
# interior(model.tracers.T, rand(1:10), rand(1:10), :), p)
201-
# @test isequal(trues(length(ρ_profile)),
202-
# test_density_profile(d_field, ρ_profile, atol))
203-
# end
196+
end

0 commit comments

Comments
 (0)