Changes to naming in the TurbulenceClosures module

[tomchor]

This PR will rename some variables and funcitons in the TurbulenceClosures module to make things more clear.

It will also remove the fallback dispatch of calc_ κᶜᶜᶜ() and define calc_ κᶜᶜᶜ() at least for the Smagorisnky-Lilly closure.

It also changes the behavior of viscosity() to return a tuple when given a tuple, instead of summing all the viscosities in the tuple. This avoids misuse by unattentive users who have a tuple with different-formulation-viscosities (for example horizontal and vertical formulations). (This comment provides an example.)

Closes #2751

CC @simone-silvestri

[simone-silvestri]

We could also just change the signature of the calc_ functions to pass all the fields instead of just velocities and tracers separately fields(model)

[tomchor]

We could also just change the signature of the calc_ functions to pass all the fields instead of just velocities and tracers separately fields(model)

Yeah I like that idea. Also, is there anything that keeps us from passing model itself? It seems like it'd simplify the interfaces

[simone-silvestri]

There is no adapt method for the model. I don't know how it would play to send the whole model to the GPU. We might have some parameter space issue

[simone-silvestri]

I am not sure this is needed, κₑ is given as a BinaryOperation so calculating νₑ should be enough

[tomchor]

My thinking here is that it'd be good to unify the interface across models. I was also planning on creating dispatches of this function for ScalarDiffusivity() and tuple closures. It'd make things like these easier (in Oceanostics): https://github.com/tomchor/Oceanostics.jl/blob/29a544d3decd833d9f86da05b66a7392197c3b93/src/Oceanostics.jl#L25-L39

For example, as a user, I don't know what function I should use to get the diffusivities and viscosities that works for different closures. As far as I know there isn't one. This is my attempt to get us closer.

Am I missing something?

[simone-silvestri]

Ok I agree, but I would also leave that calc_nonlinear_κᶜᶜᶜ is used to calculate the diffusivity while κᶜᶜᶜ is used to retrieve the diffusivity

[glwagner]

Unity may be a good goal, but we shouldn't add code until the moment we need it / without tests.

[tomchor]

I was thinking we could change viscosity() and diffusivity() to account for the different (vertical, horizontal and 3D) formulations (with tests for it, ofc). In that way viscosity() would return, for example, [v, v, 0] for a horizontal formulation, [0, 0, v] for a vertical formulation, and a Number v for a 3D formulation. I think this would be desirable from the user's perspective and it would make diagnostics easier when non-3D formulations are used.

Based on @glwagner's comments I was also going to remove the extra definition of calc_κccc() for Smagorinsky that's currently in this PR.

However, I noticed this breaks some internal functions. Before I go ahead and fix those, I wanted to check if this changed is desirable from the developers' perspective. @glwagner @simone-silvestri what do you guys think?

[glwagner]

Why introduce a vector? What does it represent? (The diagonal elements of a hypothetical viscosity tensor, which is defined implicitly?) How would it be used?

[tomchor]

Why introduce a vector? What does it represent? (The diagonal elements of a hypothetical viscosity tensor, which is defined implicitly?) How would it be used?

In cases where the viscosity is a tuple, it would return the correct viscosities/diffusivities in each direction. For example in a case like the snippet below

tri = ScalarDiffusivity(ThreeDimensionalFormulation(), ν=1e-3, κ=1e-4)
hor = ScalarDiffusivity(HorizontalFormulation(), ν=1e-3, κ=(b=1e-4, c=1e-5))
ver = ScalarDiffusivity(VerticalFormulation(), ν=1e-4, κ=(b=1e-5, c=1e-6))
tup = (hor, ver, tri)

using viscosity(model.closure, model.diffusivity_fields) in this PR returns the expected viscosity vector with different viscosities in each direction. This would, imo, make calculating diagnostics where the viscosity is important (most notably dissipation terms) easier for the user.

Using the same code in the master branch returns one number (1e-3+1e-3+1e-4) as the viscosity, which isn't what one would expect, and which requires some extra attention and code from the user when calculating viscous diagnostics.

[simone-silvestri]

I see. It would be nice to have this as a user API convenience, but then to be complete we must include all the other entries of the tensor. It might get a bit tedious when you want to include isopycnal diffusivities like TwoDimensionalLeith or GM (IsopycnalSkewSymmetricDiffusivity) for which we do not compute the tensor but directly the flux

[tomchor]

I see. It would be nice to have this as a user API convenience, but then to be complete we must include all the other entries of the tensor. It might get a bit tedious when you want to include isopycnal diffusivities like TwoDimensionalLeith or GM (IsopycnalSkewSymmetricDiffusivity) for which we do not compute the tensor but directly the flux

I think the other tensor entries are much less frequently needed, so I'd propose we start only with the simple stuff and add the rest when needed.

I was also thinking that for inner-working purposes it's best to keep the current viscosity() behavior (always returning a 'Number`) and create a second function that can be exported to the user which has the array behavior I described. Thoughts?

[glwagner]

Why is a vector that represents the diagonal elements of a hypothetical viscosity tensor useful?

[tomchor]

Why is a vector that represents the diagonal elements of a hypothetical viscosity tensor useful?

Because I think most of the tuple closures used are (HorizontalDiffusivity, VerticalDiffusivity), where the diagonal is what you need. Although that intuition might be wrong.

But most of the reason for my attempted changes to viscosity() here is that apparently viscosity() is the user interface to retrieve the viscosities regardless of closure. However, if a user uses that in the example below, the output isn't correct considering the physics:

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> closure = (HorizontalScalarDiffusivity(ν=1), VerticalScalarDiffusivity(ν=2));

julia> model = NonhydrostaticModel(grid=grid, closure=closure);

julia> using Oceananigans.TurbulenceClosures: viscosity

julia> viscosity(model.closure, model.diffusivity_fields)
3.0

Maybe the best way to move forward isn't to change viscosity(), but IMO a user-facing function to get viscosities that works as expected (i.e., returns something like [1, 1, 2] in the above example) would be helpful.

[glwagner]

I think what I'm asking is what you envision using this vector for. I don't see the purpose of it. One could just as easily inspect each closure in the tuple individually to see what the viscosity for each closure type is. This method is unambiguous and more general than developing a "vector abstraction" for the diagonal components of a viscosity tensor.

That said, I agree that the output of viscosity for a closure tuple isn't useful. It's summing the viscosities together. This makes sense when the formulations are the same for each component, but not when they are different. So we should change that.

The main use case envisioned is when you have two ScalarDiffusivity with ThreeDimensionalFormulation. In that case, viscosity returns an object (a BinaryOperation) representing the sum of a the nonlinear diffusivity (a field) and the background molecular component (a number). This can be used in subsequent AbstractOperations.

When the closure tuple involves multiple ScalarDiffusivity with heterogeneous formulations, we need a different abstraction.

We also need to deal with the case where the closure tuple contains non-scalar-diffusivity closures.

To design an abstraction, I think we should start with a use case, which can help us develop requirements for the abstraction. Once we have requirements, we can implement something minimal and simple that's easy to reason with and that will generalize to more complicated use cases we may want to consider in the future (hopefully).

[tomchor]

I think what I'm asking is what you envision using this vector for. I don't see the purpose of it. One could just as easily inspect each closure in the tuple individually to see what the viscosity for each closure type is. This method is unambiguous and far more general than developing a "vector abstraction" for the diagonal components of a viscosity tensor.

That said, I agree that the output of viscosity for a closure tuple isn't useful. It's summing the viscosities together. This makes sense when the formulations are the same for each component, but not when they are different. So we should change that.

The main use case envisioned is when you have two ScalarDiffusivity with ThreeDimensionalFormulation. In that case, viscosity returns an object (a BinaryOperation) representing the sum of a the nonlinear diffusivity (a field) and the background molecular component (a number). This can then be used in subsequent AbstractOperations for calculations.

When the closure tuple involves multiple ScalarDiffusivity with heterogeneous formulations, we need a different abstraction.

We also need to deal with the case where the closure tuple contains non-scalar-diffusivity closures.

To design an abstraction, I think we should start with a use case, which can help us develop requirements for the abstraction. Once we have requirements, we can implement something minimal and simple that's easy to reason with and that will generalize to more complicated use cases we may want to consider in the future (hopefully).

I agree with this. I'll revert viscosity() to its original formulation for now. My only question/suggestion is: should we check in viscosity(::Tuple) that the all the elements have the same formulation? That way we can throw a warning (or maybe even an error) when trying to add different formulation closures.

[glwagner]

I wonder if we should return a tuple instead, and let the user compute the sum if that's valid? Eg

viscosity(diffusivities, closure_tuple) # returns a tuple

# For users who know this is valid:
nu = sum(viscosity(diffusivities, closure_tuple))

[tomchor]

I wonder if we should return a tuple instead, and let the user compute the sum if that's valid? Eg

viscosity(diffusivities, closure_tuple) # returns a tuple

# For users who know this is valid:
nu = sum(viscosity(diffusivities, closure_tuple))

I like this idea! It's implemented right now.

@glwagner @simone-silvestri Let me know if there's anything left to adjust in this PR

[glwagner]

Why this change?

[tomchor]

The rationale was to differentiate from calc_νᶜᶜᶜ() from νᶜᶜᶜ(), which at first seem to do the same thing. That said, I'm okay with other names. What would you suggest?

[glwagner]

Do we even need separate calc_νᶜᶜᶜ() and νᶜᶜᶜ() if they do the same thing?

[tomchor]

I don't know. I always assumed there was a reason to separate them. If we can use only one of them and get rid of the other, I'm all for it.

[glwagner]

I think some of the changes here are positive. However, I don't see the point of changing calc_νᶜᶜᶜ to calc_nonlinear_νᶜᶜᶜ (it's ambiguous what "nonlinear" refers to, so I don't think the name change improves the situation --- that's not to say there isn't room for improvement). I also prefer c to tracer.

[tomchor]

I think some of the changes here are positive. However, I don't see the point of changing calc_νᶜᶜᶜ to calc_nonlinear_νᶜᶜᶜ (it's ambiguous what "nonlinear" refers to, so I don't think the name change improves the situation --- that's not to say there isn't room for improvement). I also prefer c to tracer.

Thanks for the review. I'm relatively agnostic about what we call things. Mostly I'd just like to make this module more readable.

[glwagner]

I think the main issue here is that there is too much code, reflecting the fact that an interface for defining closures has emerged over time rather than being designed from the ground up. It could possibly benefit from a rethink. It's not as much an issue of the names of things in my opinion.

[tomchor]

@glwagner this almost became stale but I think it's ready to review and possibly merge. It doesn't re-work the code like you suggest here, but it does make the functions more easily understandable on first pass by being more verbose.

It also changes the behavior of viscosity() and diffusivity(), which no longer sum over all individual closures by default when the closure is a tuple, avoiding user error (as you suggested here).

Finally it also changes the function calc_κᶜᶜᶜ ro calc_nonlinear_κᶜᶜᶜ in order to differentiate from κᶜᶜᶜ. Although I'm agnostic about the name and could change it to anything else.

This PR used to also remove a fallback that was a bit problematic because it silently returned zero diffusivity values for closures like Smagorinsky, but it seems some other PR got around to that first before merging this one :)

[glwagner]

Ok, I will review! Except, tests are failing? Also should we merge main?

[tomchor]

Ok, I will review! Except, tests are failing? Also should we merge main?

Thanks!

And yeah there was a typo in my last commit but it should be fixed as tests should be passing. I also merged main in commit 479056a (yesterday) so we should be good regarding that.

[tomchor]

@glwagner Tests are passing now!

[glwagner]

Looks good!