Setting different diffusivities for salinity and temperature tracers

[jbisits]

The mwe (probably could be more minimal) below shows the same model/simulation setup except for whether the salinity or temperature diffusivity is set first in the closure. Given the closure is set with a NamedTuple I would have expected the same result in both cases. Instead, it seems like it is done in some kind of sequential/alphabetical order so that in the first simulation, where salinity diffusivity appears first (κ=(S=1e-7, T=1e-5)) in the closure, the salinity is the slower diffusing tracer (as expected)

but in the second simulation where the temperature diffusivity is set first (κ=(T=1e-5, S=1e-7)) the salinity is the faster diffusing tracer even though it should still be the slower diffusing tracer

Is this the expected behaviour? No problem if not but I could not see anything in the docstring about specifying the order of the diffusivities so might be worth adding something? Or is there something else at play here?

MWE:

using Oceananigans, CairoMakie

## Setting S diffusivity first
grid = RectilinearGrid(CPU(),
                       topology = (Periodic, Flat, Bounded),
                       size = (5, 50),
                       x = (-0.01/2, 0.01/2),
                       z = (-1, 0))

buoyancy = SeawaterBuoyancy()
tracers = (:S, :T)

closure = ScalarDiffusivity(; ν=1e-4, κ=(S=1e-7, T=1e-5)) # S diff first

model = NonhydrostaticModel(; grid, buoyancy, tracers, closure)
Tᵤ, Tₗ = -1.5, 0.5
T₀(x, z) = z <= -0.5 ? Tᵤ : Tₗ
Sᵤ, Sₗ = 34.58, 34.57
S₀(x, z) = z <= -0.5 ? Sᵤ : Sₗ
set!(model, S = S₀, T = T₀)

simulation = Simulation(model, Δt=0.1, stop_time=60*60)
wizard = TimeStepWizard()
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

tracers = Dict("S" => model.tracers.S, "T" => model.tracers.T)
simulation.output_writers[:tracers] = JLD2Writer(model, tracers;
                                                filename = "tracers_κₛ_κₜ",
                                                schedule = TimeInterval(60))

run!(simulation)

S = FieldTimeSeries("tracers_κₛ_κₜ.jld2", "S")
T = FieldTimeSeries("tracers_κₛ_κₜ.jld2", "T")
x = xnodes(S.grid, Center())
z = znodes(S.grid, Center())

fig = Figure(size = (1000, 1000))
axS₀ = Axis(fig[1, 1], ylabel = "z", title = "Salinity t = 0")
heatmap!(axS₀, x, z, interior(S, :, 1, :, 1), colormap = :haline)
axT₀ = Axis(fig[1, 2], ylabel = "z", title = "Temperature t = 0")
heatmap!(axT₀, x, z, interior(T, :, 1, :, 1), colormap = :thermal)
axSend = Axis(fig[2, 1], xlabel = "x", ylabel = "z", title = "Salinity t = 61")
heatmap!(axSend, x, z, interior(S, :, 1, :, 61), colormap = :haline)
axTend = Axis(fig[2, 2], xlabel = "x", ylabel = "z", title = "Temperature t = 61")
heatmap!(axTend, x, z, interior(T, :, 1, :, 61), colormap = :thermal)
fig
save("Sdiff_set_first.png", fig)

## Setting T diffusivity first
grid = RectilinearGrid(CPU(),
                       topology = (Periodic, Flat, Bounded),
                       size = (5, 50),
                       x = (-0.01/2, 0.01/2),
                       z = (-1, 0))

buoyancy = SeawaterBuoyancy()
tracers = (:S, :T)

closure = ScalarDiffusivity(; ν=1e-4, κ=(T=1e-5, S=1e-7)) # T diff set first

model = NonhydrostaticModel(; grid, buoyancy, tracers, closure)
Tᵤ, Tₗ = -1.5, 0.5
T₀(x, z) = z <= -0.5 ? Tᵤ : Tₗ
Sᵤ, Sₗ = 34.58, 34.57
S₀(x, z) = z <= -0.5 ? Sᵤ : Sₗ
set!(model, S = S₀, T = T₀)

simulation = Simulation(model, Δt=0.1, stop_time=60*60)
wizard = TimeStepWizard()
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

tracers = Dict("S" => model.tracers.S, "T" => model.tracers.T)
simulation.output_writers[:tracers] = JLD2Writer(model, tracers;
                                                filename = "tracers_κₜ_κₛ",
                                                schedule = TimeInterval(60))

run!(simulation)

S = FieldTimeSeries("tracers_κₜ_κₛ.jld2", "S")
T = FieldTimeSeries("tracers_κₜ_κₛ.jld2", "T")
x = xnodes(S.grid, Center())
z = znodes(S.grid, Center())

fig = Figure(size = (1000, 1000))
axS₀ = Axis(fig[1, 1], ylabel = "z", title = "Salinity t = 0")
heatmap!(axS₀, x, z, interior(S, :, 1, :, 1), colormap = :haline)
axT₀ = Axis(fig[1, 2], ylabel = "z", title = "Temperature t = 0")
heatmap!(axT₀, x, z, interior(T, :, 1, :, 1), colormap = :thermal)
axSend = Axis(fig[2, 1], xlabel = "x", ylabel = "z", title = "Salinity t = 61")
heatmap!(axSend, x, z, interior(S, :, 1, :, 61), colormap = :haline)
axTend = Axis(fig[2, 2], xlabel = "x", ylabel = "z", title = "Temperature t = 61")
heatmap!(axTend, x, z, interior(T, :, 1, :, 61), colormap = :thermal)
fig
save("Tdiff_set_first.png", fig)

[glwagner]

this is a bug for sure

[glwagner]

this should be quick, but it might be helpful if you can add a test. I will open a PR soon

[glwagner]

actually I am not sure that causes the issue yet --- standby

[glwagner]

hmm ok the plot thickens. Here is an attempt at an MWE which shows no issue:

using Oceananigans
grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
closure = ScalarDiffusivity(κ=(a=1, b=2))
model = NonhydrostaticModel(; grid, closure, tracers=(:a, :b))
@show model.closure

model = NonhydrostaticModel(; grid, closure, tracers=(:b, :a))
@show model.closure

The problem arises later I think. There is an assumption that the named tuple for diffusivities has the same ordering as the tracers.

[glwagner]

ok found it, fixed here: #4437

[glwagner]

should be fixed*

[jbisits]

Thanks @glwagner!!

[jbisits]

I can close either close this now or wait till #4437 is merged to do so --- let me know what you think is best.

[navidcy]

Given that it's a discussion it can't be closed; only mark what provided the solution (which already has been done!)

Thanks @jbisits!