Sign error for mean PV gradient Qy in MultiLayerQG module

[mjclobo]

This is an issue made to accompany my pull request (#329 ) for this bug.

There is a sign error in the multilayerqg.jl calculation of Qy:

GeophysicalFlows.jl/src/multilayerqg.jl

Line 367 in c0c42ce

CUDA.@allowscalar @views Qy[:, :, j] = @. Qy[:, :, j] - Fp[j] * (U[:, :, j+1] - U[:, :, j]) + Fm[j-1] * (U[:, :, j-1] - U[:, :, j])

The first minus sign should be a plus sign. This creates erroneously large interior PV gradients.

The following figures can be reproduced using this script and my forked version of GeophysicalFlows.jl. For each figure the left plot is of the meridional PV gradient and the right plot is of the vertical structure of the most unstable mode. The reduced gravity error mentioned here is brought up in discussion #325 .

The first figure shows results when using density values of O(1), as in the Phillips model example. The sign error leads to a large erroneous positive interior PV gradient and a deep Charney-type instability that arises from interactions between the negative bottom PV gradient and the positive interior PV gradient. The reduced gravity error creates an erroneous negative interior PV gradient, that then interacts with the positive surface PV gradient to manifest a surface-intensified Charney-like instability. Fixing both errors leads to no instability since the Eady edge waves are too far apart (given the current background conditions) to interact and create an Eady-type instability.

@Wendazhang33, my office mate, pointed out that when you use realistic ocean density values, of O(1000), then the reduced gravity definition ceases to matter. In this case the orange line falls directly under the green line. Ironically, when O(rho)>>O(delta rho), the Boussinesq assumption applies, which means that a constant reference density should still be used in the denominator of the definition of reduced gravity. Regardless, the sign error is of the same order of consequence as when using O(1) density values, as is shown in the second plot, below. Also note that with realistic ocean density values, the corrected code leads to a nice Eady instability, as expected.

[navidcy]

thanks @mjclobo!!

Note, if you write your post like below (with an empty line between txt and link to code), the code will appear in the post.

this is some piece of code I'd like to attach

https://github.com/FourierFlows/GeophysicalFlows.jl/blob/c0c42cebd076b9cf828af1b956a08241a70febcc/src/multilayerqg.jl#L367-L367

I think this and that....

(You can edit your post and see.)

[navidcy]

Regarding the sign error:

Here's the code that computes Qy for each layer:

GeophysicalFlows.jl/src/multilayerqg.jl

Lines 365 to 369 in c0c42ce

	CUDA.@allowscalar @views Qy[:, :, 1] = @. Qy[:, :, 1] - Fp[1] * (U[:, :, 2] - U[:, :, 1])
	for j = 2:nlayers-1
	CUDA.@allowscalar @views Qy[:, :, j] = @. Qy[:, :, j] - Fp[j] * (U[:, :, j+1] - U[:, :, j]) + Fm[j-1] * (U[:, :, j-1] - U[:, :, j])
	end
	CUDA.@allowscalar @views Qy[:, :, nlayers] = @. Qy[:, :, nlayers] - Fm[nlayers-1] * (U[:, :, nlayers-1] - U[:, :, nlayers])

In #327 you suggested changing

- Fp[j] * (U[:, :, j+1] - U[:, :, j]) -> + Fp[j] * (U[:, :, j+1] - U[:, :, j]) for j = 2,...,n-1

But then this is inconsistent with layers 1 and n.

Do you agree with the equations in the docs. If yes, then it seems to me that the sign error is on the Fm term in the interior, i.e.,

+ Fm[j-1] * (U[:, :, j-1] - U[:, :, j] -> - Fm[j-1] * (U[:, :, j-1] - U[:, :, j]

Right?

[navidcy]

@Wendazhang33, my office mate, pointed out that when you use realistic ocean density values, of O(1000), then the reduced gravity definition ceases to matter. In this case the orange line falls directly under the green line. Ironically, when O(rho)>>O(delta rho), the Boussinesq assumption applies, which means that a constant reference density should still be used in the denominator of the definition of reduced gravity. Regardless, the sign error is of the same order of consequence as when using O(1) density values, as is shown in the second plot, below. Also note that with realistic ocean density values, the corrected code leads to a nice Eady instability, as expected.

Yes, I was realizing last night that we might be using the Boussinesq approximation under the hood in the buoyancy frequency definition. I'm trying to clear this up.

[mjclobo]

Regarding the sign error:

Here's the code that computes Qy for each layer:

GeophysicalFlows.jl/src/multilayerqg.jl

Lines 365 to 369 in c0c42ce

	CUDA.@allowscalar @views Qy[:, :, 1] = @. Qy[:, :, 1] - Fp[1] * (U[:, :, 2] - U[:, :, 1])
	for j = 2:nlayers-1
	CUDA.@allowscalar @views Qy[:, :, j] = @. Qy[:, :, j] - Fp[j] * (U[:, :, j+1] - U[:, :, j]) + Fm[j-1] * (U[:, :, j-1] - U[:, :, j])
	end
	CUDA.@allowscalar @views Qy[:, :, nlayers] = @. Qy[:, :, nlayers] - Fm[nlayers-1] * (U[:, :, nlayers-1] - U[:, :, nlayers])

In #327 you suggested changing

- Fp[j] * (U[:, :, j+1] - U[:, :, j]) -> + Fp[j] * (U[:, :, j+1] - U[:, :, j]) for j = 2,...,n-1

But then this is inconsistent with layers 1 and n.

Do you agree with the equations in the docs. If yes, then it seems to me that the sign error is on the Fm term in the interior, i.e.,

+ Fm[j-1] * (U[:, :, j-1] - U[:, :, j] -> - Fm[j-1] * (U[:, :, j-1] - U[:, :, j]

Right?

Shoot, you're right. Thanks for the catch! The results are the same for the constant shear/stratification case, but will definitely not be the same generally.

I closed the pull request. Is this the correct thing to do? Should I open a new one, or will you do it? Thanks for the help!

[navidcy]

Open a new one only with the sign error fix?

Let's leave the reference density discussion for a different PR.

[navidcy]

Try to keep the PR clean (i.e., don't include example scripts of your own in the PR as a suggestion to be merged in the main branch). You can always post code as part of a comment on the PR or here! :)

[mjclobo]

Thanks for the advice @navidcy. I just opened PR #329.

Is it okay to have the interim commits in the pull request, since the actual difference between the codes ends up only being the one line? Or should I make sure to only have commits relevant to the issue shown? I couldn't figure out how to do the latter, but am happy to keep trying if it's preferable.

[navidcy]

That's alright :) When we merge we can do "Squash and merge" and that will compress all the commits into one.

[navidcy]

PR #329 dealt with the sign error!

[navidcy]

@mjclobo note that I also wanted to revisit the reference density issue. That's partly described in #325 but this issue also includes some relevant plots.

I think it's good to have this issue closed since #329 actually dealt with the issue title.

[mjclobo]

@mjclobo note that I also wanted to revisit the reference density issue. That's partly described in #325 but this issue also includes some relevant plots.

I think it's good to have this issue closed since #329 actually dealt with the issue title.

Yes, my bad! I realized that I conflated the two separate issues here..I will add plots related to the reduced gravity definition only to #325 later today.