Option to specify PV gradients explicitly rather than computing them spectrally

[apaloczy]

In a two-layer MultiLayerQG configuration with topography based on the Phillips problem example from the docs, specifying the topographic PV gradient analytically instead of computing it spectrally avoids numerical problems at the boundaries for a planar bottom slope. This seems to be because of the discontinuities introduced by the FFTs.

So an option to specify PV gradients like etax and etay directly for simple idealized topographies could be helpful. A while ago I also experimented with specifying Uyy analytically for a Bickley jet in the upper layer. Maybe that's another desirable feature that could be implemented in the same way.

Any thoughts on how to go about this? Perhaps using multiple dispatch with the Params struct? I'd be happy to open a PR once we decide how to implement it.

Example 1) Spectrally calculated $\eta_y = \mathrm{IFFT}[ il \times \mathrm{FFT}(h) ] \times f_0/H_2$ with a linear slope $h = sy$:

https://user-images.githubusercontent.com/4955404/192116919-2324ce79-67e9-4eda-9f16-0f799006daab.mp4

Example 2) Directly prescribed $\eta_y = s \times f_0/H_2$:

https://user-images.githubusercontent.com/4955404/192117079-bd1c1ce2-536e-40e2-b712-5f223c2a395f.mp4

[navidcy]

Good point! Same issue would be applicable to SingleLayerQG module. Since always is the PV gradients that come into play, we. can prescribe the topographic PV gradient as:

η_x * x + η_y * x + η(x, y)

with η_x, η_y constants denoting the large-scale topographic gradients and η a function periodic in x, y. How does that sound?

[apaloczy]

That sounds good and it would fix the linear slope cases, but would leave out other idealized topographies such as h = cos(x)cos(y). Although I suppose the error in computing the gradient spectrally should be close to machine precision in most periodic h cases?

Would having the option to prescribe either the full PV gradients η_x(x, y), η_y(x, y) or the topographic PV η(x, y) be a good approach? In the latter case we would compute the PV gradients spectrally (as is currently done), and in the former we would just add them directly to Qx and Qy.

[navidcy]

Yes, if η is periodic then the slope is correct up to eps when computed with ffts. The problem arises when you try to apply fft to a non-periodic function.

We can try it out perhaps and see if we get differences?

[apaloczy]

Right, I checked the maximum absolute differences for the derivatives of h = cos(x)cos(y) with N = 64, 128, 256, 512, 1024, and they are O(1e-14) or O(1e-13).