NonhydrostaticModel creates an extra empty element in diffusivities list when a tuple is passed to as closure

[tomchor]

I'm assuming the behavior below isn't expected., where NonhydrostaticModel creates an element in diffusivities for each element in a tuple of closures but one of them is empty. What would be the preferred way?

julia> using Oceananigans

julia> grid = RegularRectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))
RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (4, 4, 4)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.25, 0.25, 0.25)

julia> closure = (AnisotropicMinimumDissipation(), IsotropicDiffusivity(ν=1e-6, κ=1.4e-7))
(AnisotropicMinimumDissipation{Float64} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.08333333333333333
                        Buoyancy modification multiplier Cb: nothing
                Background diffusivit(ies) for tracer(s), κ: 0.0
             Background kinematic viscosity for momentum, ν: 0.0, IsotropicDiffusivity: ν=1.0e-6, κ=1.4e-7)

julia> model = NonhydrostaticModel(grid=grid, closure=closure)
NonhydrostaticModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=4, Ny=4, Nz=4)
├── tracers: (:T, :S)
├── closure: Tuple{AnisotropicMinimumDissipation{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Float64, NamedTuple{(:T, :S), Tuple{Float64, Float64}}, Float64, NamedTuple{(:T, :S), Tuple{Float64, Float64}}, Nothing}, IsotropicDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Float64, NamedTuple{(:T, :S), Tuple{Float64, Float64}}}}
├── buoyancy: SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}
└── coriolis: Nothing

julia> model.diffusivities[1]
(νₑ = Field located at (Center, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (4, 4, 4)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=4, Ny=4, Nz=4)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux, κₑ = (T = Field located at (Center, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (4, 4, 4)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=4, Ny=4, Nz=4)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux, S = Field located at (Center, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (4, 4, 4)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=4, Ny=4, Nz=4)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux))

julia> model.diffusivities[2]

Note that the output of the last line is empty. So the second element is there, but it's nothing. If I check the model's property closure, I can see both closures though:

julia> model.closure[1]
AnisotropicMinimumDissipation{Float64} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: (T = 0.08333333333333333, S = 0.08333333333333333)
                        Buoyancy modification multiplier Cb: nothing
                Background diffusivit(ies) for tracer(s), κ: (T = 0.0, S = 0.0)
             Background kinematic viscosity for momentum, ν: 0.0

julia> model.closure[2]
IsotropicDiffusivity: ν=1.0e-6, κ=(T = 1.4e-7, S = 1.4e-7)

So this is a bit confusing (at least to me) and I guess it relates to #1381.

Thoughts?

[glwagner]

For tupled closures we use a tuple of diffusivities. So the in your example, diffusivities[2] corresponds to what you'd get if you build a model with just IsotropicDiffusivity --- isnothing(model.diffusivities).

Tuple closures are certainly a bit rough... what do you expect to have here? What would you expect to happen with, for example, three closures, two of which have diffusivities but one of which does not?

[glwagner]

Part of the problem we are trying to solve is how to evaluate the closure term for multiple closure models. For this we basically unroll a loop (using recursion) over all the closures, and add each term separately. Since the closure term (the flux divergence) expects a diffusivity, we need diffusivities for each closure in a tuple.

For example, see how the closure term is evaluated for a tuple of two closures:

Oceananigans.jl/src/TurbulenceClosures/closure_tuples.jl

Lines 14 to 16 in 201b50e

	@inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{C1, C2}, clock, U, Ks, args...) where {C1, C2} = (
	$stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)
	+ $stress_div(i, j, k, grid, closures[2], clock, U, Ks[2], args...))

[tomchor]

what do you expect to have here? What would you expect to happen with, for example, three closures, two of which have diffusivities but one of which does not?

I don't really understand the question though. What do you mean by a closure without diffusivity? The way I see it every closure has a diffusivity associated with it (and here I'm counting viscosity as a momentum diffusivity), whether it's a constant diffusivity or a dynamically-calculated (eddy) one. (The exception being closure=nothing, but that's the absence of closure.)

That said, what I expected was model.diffusivities[2] to store the constant values ν=1e-6 and κ=1.4e-7, since that's the second closure of the tuple.

[tomchor]

One other option that might be easier to implement is to keep the store the tuple in model.closures, but only store the total viscosity (the sum of all closures in the tuple) in model.diffusivities. Maybe we'd have to rename it as model.diffusivities.ν_tot and model.diffusivities.κ_tot .

I think that's as clear, but simpler.

Scratch that! I was thinking of a specific example I had in mind and forgot that not every closure can be simplified to the sum of diffusivities times the nabla operator.

[glwagner]

Ah I see the confusion. The field model.diffusivities stores viscosities and diffusivities that have to be precomputed. For diffusivities prescribed as constants or functions we don't need to store the diffusivities; instead the information needed to compute the diffusivity or viscosity is bound to model.closure.

Perhaps we could call this model property "diffusivity_fields"? This would distinguish the information there from constant or function diffusivities.

[tomchor]

Yeah I think this change would help! It would definitely make things clearer. Here's a couple alternative thoughts though: Would it be possible to get rid of diffusivities altogether and store the precomputed fields in the closures directly? E.g. AMD would have an extra field data, so would smagorinsky and so on. I think it'd simplify things, although I'm not sure how feasible that would be. I guess in order to calculate those, the closure would need to already be associated with a model, so that's a downside right away... Another thought: why not just also store the constant viscosities in the diffusivities property? I know it's not necessary, but it wouldn't take up any significant space and it would make things more intuitive/organized I think.

…

On Tue, Jul 20, 2021, 12:40 Gregory L. Wagner ***@***.***> wrote: Ah I see the confusion. The field model.diffusivities stores viscosities and diffusivities that have to be precomputed. For diffusivities prescribed as constants or functions we don't need to store the diffusivities; instead the information needed to compute the diffusivity or viscosity is bound to model.closure. Perhaps we could call this model property "diffusivity_fields"? This would distinguish the information there from constant or function diffusivities. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1878 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADEX5KSF44HXRQARFRTAQ6DTYXGJNANCNFSM5AV5LNJQ> .

[glwagner]

An additional thought: it might be nice to great an auxiliary_fields object for users and closures alike and put diffusivities there. The issue is closure tuples. For example, god forbid a user wants two LES closures. This might require two diffusivity fields with the same name (eg nu_e). Then this method fails. We could have all closures use unique names, like AnisotropicMinimimDissipation_nu_e. But then I'm wondering if this is much of an improvement to the code, since it adds complexity of its own (and we have plenty of work to do).

[glwagner]

Those are all possible developments.

To put fields in closure we have to implement a pattern (implemented for things like boundary conditions, though it's taken a lot of work to get it right) whereby users instantiate a "template" closure, to which fields / data are added later given knowledge of the grid in the model constructor. Otherwise, users have to provide grid as an argument when constructing both a closure and a model. Either pattern is possible --- but implementing a "template + materialization" design is complicated, whereas requiring users to provide grid when constructing a closure changes the API.

When a complicated algorithm endows the code with more features, I feel it could be justified if the feature is good enough. But if its only purpose is to rearrange where array references are bound, it feels less worth the trade off of code maintenance and the work to implement it. Another solution is documentation.

Storing constant values in diffusivities is possible, but could be interpreted as boilerplate. The coefficients need to be stored in closure, because that's where users specify them. So copying the constants or functions into diffusivities (this would have to be done individually for all closures) requires additional code.

I think of all the changes to diffusivities, the one that makes the most sense is to eliminate that field and change the API so that grid is a required argument for closures that have fields associated with them. This is a substantial internal refactor, but perhaps someone wants to take it up?

Note: I don't this will work for memory parallel models until #1825 is resolved.

[tomchor]

Storing constant values in diffusivities is possible, but could be interpreted as boilerplate. The coefficients need to be stored in closure, because that's where users specify them. So copying the constants or functions into diffusivities (this would have to be done individually for all closures) requires additional code.

We don't have that many closures that would need this though, right? I think only IsotropicConstant, AnisotropicConstant and the biharmonic one that I see in the docs. (I think you have a TKE-based one that isn't listed in the docs also?)

I think of all the changes to diffusivities, the one that makes the most sense is to eliminate that field and change the API so that grid is a required argument for closures that have fields associated with them. This is a substantial internal refactor, but perhaps someone wants to take it up?

I like this option, but keeping the user from creating a closure without a grid seems a bit counter-productive... It's still not clear to me which one I'd vote for.

For the time being, should we do something for the problem title refers to?

[glwagner]

I think all the closures we have contain diffusivity data (constants or functions or fields) that are not in diffusivities.

The list is here:

https://github.com/CliMA/Oceananigans.jl/tree/master/src/TurbulenceClosures/turbulence_closure_implementations

Only a subset of the code is documented, so better to look to the source when considering refactoring I think.

[glwagner]

For the time being, should we do something for the problem title refers to?

What can we do that will preserve the functionality of our closure tuples? For example, we need this code to work:

Oceananigans.jl/src/TurbulenceClosures/closure_tuples.jl

Lines 14 to 16 in 201b50e

	@inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{C1, C2}, clock, U, Ks, args...) where {C1, C2} = (
	$stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)
	+ $stress_div(i, j, k, grid, closures[2], clock, U, Ks[2], args...))

[tomchor]

I guess the only non-invasive option would be to leave .closure alone and copy constant diffusivities to diffusivities. Or I guess we can do nothing and live with this since it's not that big a deal for the time being haha

[tomchor]

Just for the record, @glwagner and I decided to go with this suggestion which is the simplest for the time being. The naming is more intuitive and (with the new name) is makes sense that an empty element is created when a closure doesn't require an auxiliary field.