[BUG] Error in commutation grows drammatically with filter scale.

[dhruvbalwada]

One property of the filter is that it should commute with gradients. This is satisfied by gcm-filters when the filter scale is small, but as the scale increases, this condition breaks very dramatically (see figure below). I am on a simple cartesian grid.

Here is a notebook where I tested this on a simulation: https://gist.github.com/dhruvbalwada/0d2a6541cf584502401e8b6bfe0b325b This is not a reproducible example, but I think it presents the case easily enough.

Here I compare the gcm-filters against a regular gaussian convolution filter from scipy.

The figure below shows the impact of two different filters on commutation and sum properties as a function of filter scale. While the choice of method does not impact the preservation of the sum (right plot), but the commutation with an x-gradient gets worse with filter scale for gcm-filters (left plot).

Note that the issue becomes quite bad even for filter scales of 80-100km in a 5km resolution run, suggesting that it would be hard to correctly filter out mesoscale eddies in mesoscale permitting simulation - an important base case for gcm-filters.

In the notebook I also plot the impact of the different filters on the spectra (below cell [13]), and it seems that for larger filters scales the gcm-filters does not filter at the grid scale. This grid scale variance (or noise) is the likely culprit for breakdown of the filter-commutation property.

@rabernat mentioned that there might be a solution to this coming through soon, and that I should highlight this in an issue.

[dhruvbalwada]

Maybe #134 will fix this - @NoraLoose ?

[iangrooms]

I think the bump in the spectrum at the small-scale end is probably a result of numerical instability; you could try promoting to double precision to fix that, and/or then also try factoring the Gaussian filter.

The default n_steps is chosen so that the max absolute error between the target filter and the approximation is about 0.01, which is what you saw when you used the semi-log plot in your notebook. If you need more accuracy you can increase n_steps. If you use really large n_steps you'll start to see numerical instability, in which case you should promote to double precision and factor the Gaussian.

The width of the kernel is 2 * n_steps + 1. You could choose n_steps to match the width of the gcm-filters kernel to the width of the scipy Gaussian convolution kernel. That would make the computational cost of the two methods similar, and might make the accuracy more similar too.

[dhruvbalwada]

Thanks for these suggestions Ian. I was able to try almost all, and things worked much better. I am not sure how to do the factoring of the Gaussian that you mentioned. Shown in notebook here: https://gist.github.com/dhruvbalwada/03630f1d63b0fa6d54410b22732b1fae

Setting the kernel width to be comparable to the truncation width of the scipy Gaussian conv kernel in gcm-filter worked even better than the scipy function in terms of accuracy of spectra and eventual errors in commutation of gradients. However, based on a rough test scipy did seem much faster (but this might be because I only tested on one snapshot and gcm-filters might have greater initial overhead) Here is a figure comparing the filters:

I set the n_steps as np.ceil(2*sigma*4), where sigma = Lfilter/dx/np.sqrt(12). This is based on scipy doc, which says that scipy truncates at a standard deviation that is 4 times the standard deviation (I hope I did this right).

[iangrooms]

OK, glad it's working better for you now. Closing this because I don't think it's a bug.

[NoraLoose]

Hi @dhruvbalwada, thanks for trying and testing GCM-Filters!

I am not sure how to do the factoring of the Gaussian that you mentioned.

The factoring of the Gaussian is described here and in this example.

Please let us know if you have other issues or questions, and I will reopen this issue.