Vector Laplacian on C-grid #34

NoraLoose · 2021-04-23T17:08:17Z

Summary

This PR adds a vector Laplacian on a C-grid mimicking the MOM6 implementation of viscosity, based on the friction operator suggested by Griffies & Hallberg, 2000.

Tests

While our scalar Laplacians are tested to conserve the area integral, the new vector Laplacian is designed to preserve angular momentum. A corollary of angular momentum conservation is that the vector Laplacian is invariant under solid body rotations. The latter is the test that I implemented for the vector Laplacian. Shoutout to @adcroft, who helped with designing this vector Laplacian test!

Numerical instability

While the viscosity-based filter (i.e., the filter using the vector Laplacian) seems to produce reasonable filtered fields (see also figure below), it unfortunately seems way more sensitive to numerical instability than the diffusion-based filter (the filter using the scalar Laplacian). Tests with NeverWorld2 data suggest that the old stage order (before #32) is more stable for the vector Laplacian. I don't know why... I am happy to revert this PR to "draft" until we have figured this out.

A few other (minor) notes:

For now, filter.py has two separate functions:_create_filter_func and _create_filter_func_vec; as well as methods: apply and apply_to_vector, to treat diffusion- vs. viscosity-based filters. There is probably a more elegant way to do this. Suggestions are welcome!
Naming convention of grid variables follows MOM6.
I implemented anisotropy aligned with the grid directions. Work left for a future PR: A few minor additions would allow anisotropy in arbitrary direction (aligned with spatially varying normal n=(n_1, n_2) that is input by the user) a la Smith and McWilliams, 2003 - see also MOM6 code and this issue on how to correctly implement this.

…s into vec-Laplacian-Cgrid

codecov-commenter · 2021-04-23T17:10:32Z

Codecov Report

Merging #34 (289a0ff) into master (34449a6) will increase coverage by 0.27%.
The diff coverage is 99.14%.

@@            Coverage Diff             @@
##           master      #34      +/-   ##
==========================================
+ Coverage   98.19%   98.47%   +0.27%     
==========================================
  Files           7        7              
  Lines         499      719     +220     
==========================================
+ Hits          490      708     +218     
- Misses          9       11       +2

Flag	Coverage Δ
unittests	`98.47% <99.14%> (+0.27%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
gcm_filters/kernels.py	`98.87% <97.01%> (-1.13%)`	⬇️
gcm_filters/filter.py	`97.15% <100.00%> (+0.63%)`	⬆️
tests/test_filter.py	`100.00% <100.00%> (ø)`
tests/test_kernels.py	`100.00% <100.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 34449a6...289a0ff. Read the comment docs.

NoraLoose · 2021-04-23T17:15:25Z

The figures above compare diffusion-filtered vs. viscosity-filtered velocity for NeverWorld2 data. The fields are filtered from 1/16 degree to 1 degree. As expected, the differences between the two filters become more visible as you move to higher latitudes. Close to the equator, the filters produce the same result, because here the NeverWorld2 grid is essentially Cartesian.

iangrooms · 2021-04-23T18:13:19Z

I discussed this with @NoraLoose yesterday and she suggested that the numerical stability could be related to the fact that the lat/lon Neverworld2 grid is extremely anisotropic near the poleward boundaries. I support merging this PR and then continuing to work on numerical stability in the future. It would be nice to see if the stability problems persist on a more realistic (e.g. Mercator) grid, with cells that are closer to unit aspect ratio.

NoraLoose · 2021-04-23T18:32:15Z

the numerical stability could be related to the fact that the lat/lon Neverworld2 grid is extremely anisotropic near the poleward boundaries

On the other hand, the diffusion-based filter also knows about the extreme anisotropy of the NeverWorld2 grid. Nonetheless, the diffusion-based filter works okay.

Another possibility is the handling of boundaries: this is where the diffusion-based and viscosity-based filter differ quite a bit. You can even see that in the figures above (right-most panels). There seems to be some noise building up, particularly at the southern and northern boundaries. One could potentially look into whether there is a better way to handle the wet masks than done here:

https://github.com/NoraLoose/gcm-filters/blob/7f9c8bc81ef17a291e6fdc091b6dc9906bcc9b67/gcm_filters/kernels.py#L342-L345

rabernat · 2021-04-30T15:15:45Z

Sorry for taking so long to respond to this. I was on vacation last week and have still not managed to catch up! 🙃

I have just created the group @ocean-eddy-cpt/gcm-filters-dev which now has admin rights on this repo. It includes myself, @NoraLoose and @iangrooms, but we can add anyone else interested in maintaining.

At this point, I don't want to be a bottleneck on progress here. As long as someone from @ocean-eddy-cpt/gcm-filters-dev reviews and approves the PR, they should feel empowered to merge.

I can try to review early next week if desired.

rabernat

Some superficial comments.

rabernat · 2021-04-30T15:16:21Z

gcm_filters/filter.py

@@ -83,6 +83,7 @@ def _compute_filter_spec(
            else:  # ndim==2
                n_steps = np.ceil(6.4 * filter_scale / dx_min).astype(int)

+    print("n_steps = %i" % n_steps)


No print statements in the code.

Thanks for catching this! The print statement is now removed.

rabernat · 2021-04-30T15:16:59Z

tests/test_filter.py

+    # check conservation under solid body rotation: u = cos(lat), v=0;
+    # this test fails due to numerical instabilities
+    # data_u = np.cos(geolat_u/360*2*np.pi)
+    # data_v = np.zeros_like(data_u)
+    # da_u = xr.DataArray(data_u, dims=["y", "x"])
+    # da_v = xr.DataArray(data_v, dims=["y", "x"])
+    # filtered_u, filtered_v = filter.apply_to_vector(da_u, da_v, dims=["y", "x"])
+    # xr.testing.assert_allclose(filtered_u, 0.0, atol=1e-12)
+    # xr.testing.assert_allclose(filtered_v, 0.0, atol=1e-12)


Why commented-out code? Probably best to avoid this.

I fixed the test and uncommented it.

(I had wrongly assumed that numerical instability had caused the test to fail, but it was actually a small bug in the test - which is now fixed.)

- correct typo in test that checks invariance of viscosity-based filter under solid body rotation - test was previously commented out, but passes now

- change field_u_bar to ufield_bar and similar

- If the user provides area_u or area_v that contain zeros (as is the case for e.g., MOM6 output over land) then the previous version did a division by zero; - This issue is now fixed via a conditional operation.

NoraLoose · 2021-05-25T17:03:46Z

I made a small modification to the code. The previous version did a division by zero if the user provides area_u, area_v that contain zeros (as is the case for e.g., MOM6 output over land). This is now fixed via a conditional statement.

self.recip_area_u = np.where(self.area_u > 0, 1 / self.area_u, 0)

The code still raises a RuntimeWarning: divide by zero encountered in true_divide if the user runs the vector Laplacian, but only on CPUs (using the above line with np.where) - not on GPUs (using the above line with cp.where). Weird, and not optimal.

The following version would avoid the warning on CPUs, but does not work on GPUs:

self.recip_area_u = np.divide(1, self.area_u, out=np.zeros_like(self.area_u), where=self.area_u!=0)

Can anyone think of a practical solution here? Maybe @jbusecke?

I also added a new test test_zero_area in test_kernels.py which reproduces the current warning (but checks that the division by zero does not actually happen in the Laplacian).

rabernat · 2021-05-25T17:07:19Z

Good catch! Ready for another review here?

I don't know what's up with the RTD build 🤔

NoraLoose · 2021-05-25T17:07:26Z

I'm also not sure why the Read the Docs build suddenly fails. In my last three commits I haven't changed anything that should affect the docs. Any ideas?

rabernat · 2021-05-25T17:11:28Z

Any ideas?

The most likely reason is the one of the packages that the docs depend on has released a new version with an incompatibility.

The quick solution is to figure out which one it is and then pin a specific version number here:

gcm-filters/docs/environment.yml

Lines 2 to 22 in 64b3b8e

    
           dependencies: 
        
             - python=3.8 
        
             - numpy 
        
             - scipy 
        
             - xarray 
        
             - matplotlib 
        
             - numpydoc 
        
             - sphinx 
        
             - ipython 
        
             - nbconvert 
        
             - nbformat 
        
             - ipykernel # not strictly necessary but this is nice to run notebooks in this env to test 
        
             - pandoc 
        
             - pip 
        
             - pip: 
        
               - jupyter_client 
        
               - sphinx-book-theme 
        
               - sphinx-copybutton 
        
               - sphinxcontrib-srclinks 
        
               - myst-parser 
        
               - myst-nb

Unfortunately the error does not immediately make it obvious which one. Myst-nb had a minor release (0.12.3) about a month ago, so that could be it: https://github.com/executablebooks/MyST-NB/releases/
You could try pinning myst_nb==0.12.2

NoraLoose · 2021-05-25T17:11:34Z

Good catch! Ready for another review here?

Sure, thanks! 👍

* Goal: Get the RTD build to work again

rabernat

This is awesome Nora. It will be so extremely useful in many different scenarios.

My main concern is that we are overloading the BaseLaplacian class a little too heavily. See detailed comments below.

rabernat · 2021-05-25T17:16:09Z

gcm_filters/kernels.py

@@ -290,6 +299,115 @@ def __call__(self, field: ArrayType):
 ALL_KERNELS[GridType.TRIPOLAR_POP_WITH_LAND] = POPTripolarLaplacianTpoint


+@dataclass
+class VectorLaplacian(BaseLaplacian):


Suggested change

class VectorLaplacian(BaseLaplacian):

class CgridVectorLaplacian(BaseLaplacian):

Maybe this name would be a bit more specific?

I also wonder whether we should use different base classes, e.g. BaseScalarLaplacian, BaseVectorLaplacian. This would clarify a few ambiguities in the code.

rabernat · 2021-05-25T17:19:25Z

gcm_filters/kernels.py

+        self.recip_area_u = np.where(self.area_u > 0, 1 / self.area_u, 0)
+        self.recip_area_v = np.where(self.area_v > 0, 1 / self.area_v, 0)
+
+    def __call__(self, ufield: ArrayType, vfield: ArrayType):


My concern here is the the signature of this Laplacian does not match the signature of BaseLaplacian, which is

def __call__(self, field):

I'm surprised that mypy did not complain about that. In any case, it supports the idea suggested above that we should have a separate base class for vector laplacians.

Good point. Definitely a good idea to have two separate base classes.

rabernat · 2021-05-25T17:21:56Z

tests/test_filter.py

@@ -67,7 +67,11 @@ def test_filter_spec(filter_args, expected_filter_spec):
    # TODO: check other properties of filter_spec?


-@pytest.fixture(scope="module", params=list(GridType))
+vector_grids = [gt for gt in GridType if gt.name.startswith("VECTOR")]


Parsing the name string does not seem like a robust way to detect vector laplacians. Instead we should somehow refactor our base data structures to store this information. In the meantime, simply hard-coding the grid types here would not be a bad way to go. There are not so many.

rabernat · 2021-05-25T17:23:49Z

gcm_filters/filter.py

+        # with regular arguments
+        assert len(args) == len(Laplacian.required_grid_args())
+        grid_vars = {k: v for k, v in zip(Laplacian.required_grid_args(), args)}
+        laplacian = Laplacian(**grid_vars)


What happens if the user passes a non-vector Laplacian here? It should not be allowed.

rabernat · 2021-05-25T17:25:02Z

gcm_filters/filter.py

+    def apply_to_vector(self, ufield, vfield, dims):
+        """Filter a vector field across the dimensions specified by dims."""
+
+        filter_func_vec = _create_filter_func_vec(self.filter_spec, self.Laplacian)


We should be able to raise an informative error if a user tries to do a vector filter with a scalar laplacian I think we can resolve this by implementing the suggestion to have a separate BaseScalarLaplacian and BaseVectorLaplacian and then checking the type in the Filter.__post_init__ method.

Here is an attempt to raise an error based on the two base classes:
https://github.com/NoraLoose/gcm-filters/blob/59bb142a8c01224138d00026000a1732873271de/gcm_filters/filter.py#L330-L334

Does this look good to you?

- one for scalar Laplacians, one for vector Laplacians

- and vice versa

rabernat

Fantastic. Really excited to try this.

rabernat · 2021-05-25T21:06:55Z

gcm_filters/filter.py

+        if not issubclass(self.Laplacian, BaseScalarLaplacian):
+            raise ValueError(
+                f"Provided Laplacian {self.Laplacian} is a vector Laplacian. "
+                f"The .apply method is only suitable for scalar Laplacians."


Suggested change

f"The .apply method is only suitable for scalar Laplacians."

f"The ``.apply`` method is only suitable for scalar Laplacians."

rabernat · 2021-05-25T21:07:06Z

gcm_filters/filter.py

+        if not issubclass(self.Laplacian, BaseVectorLaplacian):
+            raise ValueError(
+                f"Provided Laplacian {self.Laplacian} is a scalar Laplacian. "
+                f"The .apply_to_vector method is only suitable for vector Laplacians."


Suggested change

f"The .apply_to_vector method is only suitable for vector Laplacians."

f"The ``.apply_to_vector`` method is only suitable for vector Laplacians."

rabernat · 2021-05-25T21:11:04Z

tests/test_kernels.py

+    vector_grid_type_field_and_extra_kwargs,
+):
+    """This test checks that vector Laplacians are invariant under solid body rotations:
+    a corollary of conserving angular momentum."""


Btw, I love this test! 🎉

rabernat · 2021-05-26T02:55:55Z

I manually resolved the merge conflict and it made pre-commit mad. But I'm just going to merge this so we can move forward.

NoraLoose added 24 commits March 4, 2021 14:02

initial vector Laplacian commits

dc0caa1

fix typo in grid info

d7c36fe

add apply_to_vector method to Filter class

2680219

fix bugs in vector laplacian

8bcd333

add kappas for anisotropic diffusion

a266bd7

add wet masks to vector Laplacian and fix bugs

ed5ee1a

add anisotropy to vector Laplacian

2a0a271

modify anisotropic kappas in vector Laplacian

1f8b44c

Re-order filter stages for scalar and vector Laplacian

a19f939

initial vector Laplacian commits

be9fb9a

fix typo in grid info

ca5d46a

add apply_to_vector method to Filter class

dca885d

fix bugs in vector laplacian

5766c4a

add kappas for anisotropic diffusion

175bf1a

add wet masks to vector Laplacian and fix bugs

4844be4

add anisotropy to vector Laplacian

9f85791

modify anisotropic kappas in vector Laplacian

2ff3d46

Re-order filter stages for scalar and vector Laplacian

b3e7790

Merge branch 'vec-Laplacian-Cgrid' of github.com:NoraLoose/gcm-filter…

1b39e4f

…s into vec-Laplacian-Cgrid

Remove duplicate vector Laplacian from kernels.py

e27cc7b

Update description of _create_filter_func_vec()

d10bddc

Clarify meaning of grid input variables for vector Laplacian

eb4eca4

Finalize vector Laplacian tests

d19021e

Finalize viscoscity-based filter tests

7f9c8bc

rabernat reviewed Apr 30, 2021

View reviewed changes

NoraLoose added 4 commits May 2, 2021 12:30

Correct test for viscosity-based filter

2d3e3a7

- correct typo in test that checks invariance of viscosity-based filter under solid body rotation - test was previously commented out, but passes now

Delete print statement

560ab86

Some cosmetic name changes in filter.py

396328f

- change field_u_bar to ufield_bar and similar

Delete irrelevant comment

89137f3

NoraLoose added the kernels Related to low-level kernels label May 19, 2021

NoraLoose added 3 commits May 24, 2021 20:39

Avoid division by zero in vector Laplacian

1475295

- If the user provides area_u or area_v that contain zeros (as is the case for e.g., MOM6 output over land) then the previous version did a division by zero; - This issue is now fixed via a conditional operation.

Change syntax to be compatible with GPU's

786b521

Add test that checks that no division by zero occurs

b74cf95

NoraLoose added 2 commits May 25, 2021 11:20

Pin previous version of myst_nb in environment file

9814b2e

* Goal: Get the RTD build to work again

Use correct syntax for pip dependencies

8ab0bc0

rabernat reviewed May 25, 2021

View reviewed changes

rabernat mentioned this pull request May 25, 2021

Obscure build error with 0.12.3 executablebooks/MyST-NB#330

Closed

NoraLoose added 4 commits May 25, 2021 12:05

Introduce two separate base classes for Laplacians

3fad7ea

- one for scalar Laplacians, one for vector Laplacians

Raise error if vector Laplacian is used for scalar filtering

e79b396

- and vice versa

Add test that checks that error is raised

effdc08

Hard-code definition of vector grids

59bb142

rabernat approved these changes May 25, 2021

View reviewed changes

NoraLoose and others added 3 commits May 25, 2021 15:25

Add quote characters to error message

843b75e

Merge branch 'master' into vec-Laplacian-Cgrid

d23e501

Merge branch 'master' into vec-Laplacian-Cgrid

289a0ff

rabernat merged commit ca47950 into ocean-eddy-cpt:master May 26, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Vector Laplacian on C-grid #34

Vector Laplacian on C-grid #34

NoraLoose commented Apr 23, 2021

codecov-commenter commented Apr 23, 2021 •

edited

Loading

NoraLoose commented Apr 23, 2021 •

edited

Loading

iangrooms commented Apr 23, 2021

NoraLoose commented Apr 23, 2021

rabernat commented Apr 30, 2021

rabernat left a comment

rabernat Apr 30, 2021

NoraLoose May 2, 2021

rabernat Apr 30, 2021

NoraLoose May 2, 2021

NoraLoose commented May 25, 2021

rabernat commented May 25, 2021

NoraLoose commented May 25, 2021

rabernat commented May 25, 2021

NoraLoose commented May 25, 2021

rabernat left a comment

rabernat May 25, 2021

rabernat May 25, 2021

NoraLoose May 25, 2021

rabernat May 25, 2021

rabernat May 25, 2021

rabernat May 25, 2021 •

edited

Loading

NoraLoose May 25, 2021

rabernat left a comment

rabernat May 25, 2021

rabernat May 25, 2021

rabernat May 25, 2021

rabernat commented May 26, 2021

	class VectorLaplacian(BaseLaplacian):
	class CgridVectorLaplacian(BaseLaplacian):

	f"The .apply method is only suitable for scalar Laplacians."
	f"The ``.apply`` method is only suitable for scalar Laplacians."

	f"The .apply_to_vector method is only suitable for vector Laplacians."
	f"The ``.apply_to_vector`` method is only suitable for vector Laplacians."

Vector Laplacian on C-grid #34

Vector Laplacian on C-grid #34

Conversation

NoraLoose commented Apr 23, 2021

Summary

Tests

Numerical instability

A few other (minor) notes:

codecov-commenter commented Apr 23, 2021 • edited Loading

Codecov Report

NoraLoose commented Apr 23, 2021 • edited Loading

iangrooms commented Apr 23, 2021

NoraLoose commented Apr 23, 2021

rabernat commented Apr 30, 2021

rabernat left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

NoraLoose commented May 25, 2021

rabernat commented May 25, 2021

NoraLoose commented May 25, 2021

rabernat commented May 25, 2021

NoraLoose commented May 25, 2021

rabernat left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rabernat May 25, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rabernat left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rabernat commented May 26, 2021

codecov-commenter commented Apr 23, 2021 •

edited

Loading

NoraLoose commented Apr 23, 2021 •

edited

Loading

rabernat May 25, 2021 •

edited

Loading