ocean-eddy-cpt / gcm-filters

Diffusion-based Spatial Filtering of Gridded Data
https://gcm-filters.readthedocs.io/
Other
37 stars 24 forks source link

Implementation of simple fixed factor filters suboptimal for non-centered data #165

Open NoraLoose opened 4 months ago

NoraLoose commented 4 months ago

The following issue was pointed out by @pperezhogin.

Issue

Our simple fixed factor filters work well for data centered around 0 (like velocity), but not well for fields like density. For the latter, the filter does not reduce min/max values near the boundaries and in regions of strong grid stretching.

Reason

Our simple fixed factor filters simply do

filtered_field = filter(field * area) / area

The reason why this operation does not work well for non-centered data is that it does not commute with adding a constant:

filter((field + C) * area) / area = filter(field * area) / area + C * filter(area) / area

If C is huge, like C = 1000, filtering density vs. density - 1000 will make a big difference.

Possible fix

Can we simply define the fixed factor filters as

filtered_field = filter(field * area) / filter(area)

The latter operation would commute with adding a constant.

iangrooms commented 4 months ago

This is a known property of the simple fixed-factor filters; see the penultimate sentence of section 2.6 in our paper. The proposed fix would make the filter preserve constant vectors, but it would not longer preserve the integral of the filtered field. In the case of filtering density, the filtered field would have a different total mass.

The proposed fix is sort of like saying that the filtered quantities don't live on the original grid, they live on a new grid whose areas are the filtered areas. If we could find such a new grid then the total mass would be conserved. But I think it would not be worth the effort to try to find such a new grid; it might not even be always possible.

Would it work for @Pperezhogin to compute the mean density, subtract it, filter the deviations, and then add the mean back in?

Pperezhogin commented 4 months ago

Regarding filter

filtered_field = filter(field * area) / filter(area),

it does not conserve the integral.

To @iangrooms, I currently avoid this filter (and use other filters) because it is not monotonic even after subtraction of the mean field.

Possible solution will be to introduce weightening directly into the Laplacian operator. Currently the building block is: div(grad(field area)) / area while moving area out from the grad will resolve the problem: div(areagrad(field)) / area

However, I admit that this will increase the number of operations and will make this filter closer to not simple fixed factor filter.

iangrooms commented 4 months ago

I hadn't thought of trying div(area * grad(field)) / area. Interesting idea; I'll have to think more about it. In the meantime, is the "not-simple" fixed factor filter working for you?