Update Smagorinsky-Lilly model

[tomchor]

Closes #1907

[tomchor]

I realized that Julia's \Upsilon is pretty much identical to a Y, so I changed the docs to match the code instead

[tomchor]

I made one extra change which is that I don't call N^2 / S^2 a Richardson number anymore. Since the denominator is a strain rate, not the vertical shear, I think it's odd to call it Ri. Also that variable (whatever you call it) only appears 3 times briefly so there's no need to give it a name.

It appears that our buoyancy modification doesn't really match Pressel's or Lilly's since we use an arbitrary C_b while they proposed to use 1/Pr_t. @glwagner Is there any reference for that?

[glwagner]

I don't mind if we call it "Ri" or not. The Richardson number has several interpretations and formulations; one interpretation is that it generally measures the ratio between potential energy and kinetic energy in the flow (this physical interpretation motivates its use in the "stability correction" proposed by Lilly 1962). The "rate of strain definition" limits to the "typical" oceanographic definition (also the quantity that arises in Kelvin-Helmholtz instability) for parallel laminar flow.

Lilly 1962 used "Richardson number" and the symbol "Ri" to refer to this quantity:

For me, it's a natural generalization of the Richarsdon number that applies in unidirectional / laminar flow, so I feel it's nice notation. But not using it is ok too.

[glwagner]

I made one extra change which is that I don't call N^2 / S^2 a Richardson number anymore. Since the denominator is a strain rate, not the vertical shear, I think it's odd to call it Ri. Also that variable (whatever you call it) only appears 3 times briefly so there's no need to give it a name.

It appears that our buoyancy modification doesn't really match Pressel's or Lilly's since we use an arbitrary C_b while they proposed to use 1/Pr_t. @glwagner Is there any reference for that?

I don't quite understand the question --- what would you like a reference for?

It looks like for Cb=0, the reference is the original Smagorinsky paper, whereas for Cb=1/Pr one uses the model proposed by Lilly (1962). A number in between represents something different from either --- are you wondering if there are papers out there that have used Cb other than 0 or 1/Pr? As to that I don't know.

We could change the default from Cb=1 to Cb=1/Pr or Cb=0... ? The default Pr is Pr=1 but I suppose one might inadvertently change Pr but leave Cb=1 which would lead to a different LES model than proposed by Lilly (1962).

[tomchor]

I don't quite understand the question --- what would you like a reference for? C_b=0 does not use the modification term whereas C_b=1 is the model proposed by Lilly, I think...

If I understand correctly, here's Lilly's and Pressel's model formulation:

And here's our formulation:

So it seems that our model has an extra degree of freedom that neither Lilly nor Pressel propose. To match their model, we'd need to set C_b = 1/Pr_t. I'm not arguing that that's a bad thing, but if we can't find a source for that then we possibly need to be more careful about how to present that to the user.

[tomchor]

I made one extra change which is that I don't call N^2 / S^2 a Richardson number anymore. Since the denominator is a strain rate, not the vertical shear, I think it's odd to call it Ri. Also that variable (whatever you call it) only appears 3 times briefly so there's no need to give it a name.
It appears that our buoyancy modification doesn't really match Pressel's or Lilly's since we use an arbitrary C_b while they proposed to use 1/Pr_t. @glwagner Is there any reference for that?

I don't quite understand the question --- what would you like a reference for?

It looks like for Cb=0, the reference is the original Smagorinsky paper, whereas for Cb=1/Pr one uses the model proposed by Lilly (1962). A number in between represents something different from either --- are you wondering if there are papers out there that have used Cb other than 0 or 1/Pr? As to that I don't know.

We could change the default from Cb=1 to Cb=1/Pr or Cb=0... ? The default Pr is Pr=1 but I suppose one might inadvertently change Pr but leave Cb=1 which would lead to a different LES model than proposed by Lilly (1962).

Apparently my explanation wasn't needed! Yes, that's what I mean.

I propose we either get rid of C_b and just use 1/Pr_t or make it clear in the docs and docstring that we offer an extra degree of freedom, but that the model as proposed has C_b=1/Pr_t. I'm not sure which one is best since I do like the idea of Oceananigans allowing for easy (customization), but also that favors mistakes by un-attentive users (such as myself haha).

[tomchor]

It seems like these tests are failing because they compare the LES models with some pre-computed solutions:

Oceananigans.jl/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl

Line 78 in 32c5c5a

solution₀, Gⁿ₀, G⁻₀ = get_fields_from_checkpoint(initial_filename)

If I followed the code correctly, the LES models are looped through here:

Oceananigans.jl/test/test_turbulence_closures.jl

Lines 4 to 8 in 32c5c5a

	for closure in closures
	@eval begin
	using Oceananigans.TurbulenceClosures: $closure
	end
	end

which means that those are always run with the default values. Since we're comparing with pre-computed solutions, it'd be good to explicitly specify every closure's parameters, no? The downside is that we won't be able to automatically loop through th closure like this and have to specify them by hand.

[glwagner]

So it seems that our model has an extra degree of freedom that neither Lilly nor Pressel propose.

Quite! Though Lilly admits that his model amounts to "little more than a scaling argument" and that both Pr_t and the model "constant" are themselves unknown functions. So I actually think there is room to interpret Lilly as proposing a whole family of models. Since time is finite only one has been tested I guess (and also the structure of the model is questionable, so people have moved on to other formulations...)

I propose we either get rid of C_b and just use 1/Pr_t or make it clear in the docs and docstring that we offer an extra degree of freedom, but that the model as proposed has C_b=1/Pr_t. I'm not sure which one is best since I do like the idea of Oceananigans allowing for easy (customization), but also that favors mistakes by un-attentive users (such as myself haha).

Does the default Cb = 1/Pr both retain flexibility and also reduce the chances of unexpected behavior / mistakes?

If we get rid of Cb then we may want to introduce an alternative way to eliminate the buoyancy correction entirely (currently achievable with Cb=0).

A further caveat is that using 1/Pr only applies if Pr is the same for all tracers and if a linear equation of state links tracers to buoyancy. With nonlinear equations of state it's unclear what to use for Cb. I believe we also allow a different Pr for every tracer. In this case the default Cb=1/Pr will throw an error (probably a good thing, since you'd have to decide what to use for Cb in this case).

[glwagner]

It seems like these tests are failing because they compare the LES models with some pre-computed solutions:

Oceananigans.jl/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl

Line 78 in 32c5c5a

solution₀, Gⁿ₀, G⁻₀ = get_fields_from_checkpoint(initial_filename)

~If I followed the code correctly, the LES models are looped through here: ~

Oceananigans.jl/test/test_turbulence_closures.jl

Lines 4 to 8 in 32c5c5a

	for closure in closures
	@eval begin
	using Oceananigans.TurbulenceClosures: $closure
	end
	end

which means that those are always run with the default values. Since we're comparing with pre-computed solutions, it'd be good to explicitly specify every closure's parameters, no? The downside is that we won't be able to automatically loop through th closure like this and have to specify them by hand.

True, we have to regenerate the test data to make this change to the default if we want to continue using the default in the regression test. As a quick fix we could change the regression test to use C=0.23.

[tomchor]

Does the default Cb = 1/Pr both retain flexibility and also reduce the chances of unexpected behavior / mistakes?

If we get rid of Cb then we may want to introduce an alternative way to eliminate the buoyancy correction entirely (currently achievable with Cb=0).

A further caveat is that using 1/Pr only applies if Pr is the same for all tracers and if a linear equation of state links tracers to buoyancy. With nonlinear equations of state it's unclear what to use for Cb. I believe we also allow a different Pr for every tracer. In this case the default Cb=1/Pr will throw an error (probably a good thing, since you'd have to decide what to use for Cb in this case).

Good points. Based on these points I'd argue for us to keep things as it is but include a comment about this in the docstring / docs.

Agreed?

[glwagner]

You can change the closure used in the regression test by changing this line:

Oceananigans.jl/test/test_regression.jl

Line 74 in 32c5c5a

for closure in (AnisotropicMinimumDissipation(ν=1.05e-6, κ=1.46e-7), ConstantSmagorinsky(ν=1.05e-6, κ=1.46e-7))

[tomchor]

You can change the closure used in the regression test by changing this line:

Oceananigans.jl/test/test_regression.jl

Line 74 in 32c5c5a

for closure in (AnisotropicMinimumDissipation(ν=1.05e-6, κ=1.46e-7), ConstantSmagorinsky(ν=1.05e-6, κ=1.46e-7))

Just did 2 min ago :)

[glwagner]

Does the default Cb = 1/Pr both retain flexibility and also reduce the chances of unexpected behavior / mistakes?
If we get rid of Cb then we may want to introduce an alternative way to eliminate the buoyancy correction entirely (currently achievable with Cb=0).
A further caveat is that using 1/Pr only applies if Pr is the same for all tracers and if a linear equation of state links tracers to buoyancy. With nonlinear equations of state it's unclear what to use for Cb. I believe we also allow a different Pr for every tracer. In this case the default Cb=1/Pr will throw an error (probably a good thing, since you'd have to decide what to use for Cb in this case).

Good points. Based on these points I'd argue for us to keep things as it is but include a comment about this in the docstring / docs.

Agreed?

Fine by me. I also think Cb=1/Pr would be a good choice and retains current defaults.