Define kernels for filtering vectors

[rabernat]

For models on curvilinear grids, u and v components get mixed together. So we will need to define kernels that effectively do viscosity (rather than diffusion).

[rabernat]

Now that we have the implementation worked out in more detail, we can think about how to generalize our approach to vector laplacians.

The main difference between the laplacians that we already have and the vector laplacians is that the vector ones input and output two fields instead of one.

We could imagine defining kernels like

class VectorLaplacian(BaseLaplacian):
    def __call__(self, field_u, self_field_v):
        # do stuff
        return u_component, v_component

We may also have to generalize the Filter class's apply method to account for the potentially variable number of arguments. Ideally the kernel should "know" about its expected number of inputs / ouputs, e.g. via a parameter like n_inputs.

[iangrooms]

I have a theory question about the vector Laplacian: What should it conserve? For the scalar Laplacian we want to conserve the integral, but for the vector Laplacian I'm guessing it should be something like angular momentum (which is a vector).

[rabernat]

Would be keen to get input from @StephenGriffies and @Hallberg-NOAA: they surely must think about this question in designing dissipation operators for GCMs. We basically just want to reproduce those as closely as possible in this package.

[StephenGriffies]

When considering friction in the velocity equation, it is simplest to formulate the operator as the divergence of a symmetric and trace-free stress tensor. Symmetry means the operator conserves angular momentum. Trace-free means the operator does not modify the mechanical pressure. When the viscosity is formulated in a simple manner, then the divergence of the stress tensor can be reduced to a Laplacian operator. But in general it does not get that simple when the viscosities are spatially dependent.

Apologies for the self-referencing, but I give a tutorial on these matters in Chapter 17 of my book.

[NoraLoose]

@ElizabethYankovsky: The vector Laplacian in Steve's tutorial (Chapter 17) reminds me of the Laplacian form that you were originally thinking to implement. Sorry if I diverted you and @arthurBarthe away from implementing the vector Laplacian. When we talked I thought we could get away with the simpler scalar Laplacian, even for velocities. I was wrong... Did you further look into the vector Laplacian? It would be great if we could get this to work as a group!

[ElizabethYankovsky]

@NoraLoose: the Laplacians that @arthurBarthe and I have been developing for the MOM5 kernel depend on the grid that the 'field' under consideration is on. Since it's a b-grid, the u- and v-velocity points are identical and the snippet that @rabernat showed above would be performing identical operations on the u-field and v-field. We define a separate Laplacian for fields that are oriented at T-cell centers (such as the vorticity). However, for a grid in which the u- and v-components are staggered (such as the C-grid MOM6) I agree that it is useful to have one vector Laplacian function that handles all components.

In our definition, we based the Laplacian on a 5-point stencil that is identical to the LAP_T operator in MOM5 but modified for velocity point quantities. Code for the LAP_T operator is here: https://github.com/mom-ocean/MOM5/blob/master/src/mom5/ocean_core/ocean_operators.F90

We avoided using the more complex 11-point stencil detailed in section 25.4 (and Figure 25.1) of: https://mom-ocean.github.io/assets/pdfs/MOM5_manual.pdf

In summary -- we do modify the Laplacians based on the field and its grid locations, but it's not done in 1 function as @rabernat suggests, but rather through the definition of as many Laplacians as there are quantities that reside on differing grid locations.

[NoraLoose]

@ElizabethYankovsky, the approach that you are describing is exactly the approach that I have (so far) been following too. I would label this as having 2 (or 3) scalar Laplacians on T-cells, U-cells (and V-cells). A scalar Laplacian acts only on one velocity component at a time, and does not need to know about the second velocity component. A vector Laplacian needs to know both velocity components, and acts on both of them simultaneously.

[NoraLoose]

We avoided using the more complex 11-point stencil detailed in section 25.4 (and Figure 25.1) of: https://mom-ocean.github.io/assets/pdfs/MOM5_manual.pdf

It looks like this 11-point stencil operator is exactly the vector Laplacian that we might want (and identical to the one in Steve's book).

[StephenGriffies]

Section 25.5 of the MOM5 elements

https://github.com/mom-ocean/MOM5/blob/master/doc/MOM5_elements.pdf

has the C-grid version of the friction operator as derived from taking the functional derivative of the dissipation functional and accounting for spherical geometry. The resulting expression is more general than you likely need, as the stress tensor is allowed to be anisotropic (though still symmetric and trace-free), following some work from LANL and NCAR in the late 1990s. The friction does simplify a good bit when assuming horizontal isotropy.

You are certainly going to earn your numerical methods Gold Star if you work through these details.

[iangrooms]

If we want to do filtering of velocity fields with filter scale tied to grid scale then we'll need the fully anisotropic version...

[ElizabethYankovsky]

@NoraLoose just one last point about the scalar vs. vector laplacian, apologies if I misunderstood what you wrote earlier. The way I understand this now is that there are two issues here -- one is whether we use the 11-point or the simplified 5-point stencil. The second is in whether the Laplacian is defined as a vector or a scalar operator. What you and I have done so far is applying the Laplacian operator on u and v components, which still gives us two components in the answer (but neglects certain terms that the full anisotropic Laplacian would account for). Turning this into a 'vector Laplacian' is just a coding question of whether we combine the components into one function, so that we input a vector field of u and v and output two filtered fields vs. what we have been doing so far. Did I understand this correctly? Then, the second question is whether we want to use the fully anisotropic version to more accurately filter the velocity fields.

[NoraLoose]

Turning this into a 'vector Laplacian' is just a coding question of whether we combine the components into one function, so that we input a vector field of u and v and output two filtered fields vs. what we have been doing so far.

I'm not sure if I am understanding your question correctly. Please object if I am not. The 11-point stencil lateral friction operator in Chapter 25.4 of "Elements of MOM" acts on the u- and v- components simultaneously. For instance, there are contributions from tensions and strains (and anisotropic contributions). Even if we choose the isotropic version, the new vector Laplacian has components that we don't have in the current scalar Laplacian, for instance contributions arising from horizontal shear strain. So it's not just a question of restructuring the current code, but really of reformulating it to a full vector Laplacian.

[ElizabethYankovsky]

Hi @NoraLoose, thanks for the clarification. I agree that the 11-point operator has contributions that the simplified scalar version that we've been using lacks. I'm happy to work together on formulating the 11-point version, we've discussed this a bit in the past with @LaureZanna and @arthurBarthe. The calculation involves computing each of the 12 subcell components of tension and strain (given by equations 25.35-25.58 in the MOM5 manual) -- it's a pretty intensive computation! Did you decide to completely move away from the simplified version (5-point Laplacian that we have used thus far)? Was that motivated by the fact that the full contributions from the strain and tension were missing or something else?

[NoraLoose]

@ElizabethYankovsky, I guess the discussion in #26 is what finally made us (at least me) move away from a scalar Laplacian for velocities, even though @rabernat had opened the vector Laplacian issue way before that. As long as we were filtering only on regional domains away from the Arctic and the tripolar seam, it almost seemed that we could get away with a scalar Laplacian for the velocity components. But even this is suboptimal because the u and v components should really get mixed together for models on curvilinear grids (also away from the tripolar seam), and the problem becomes worse if the filter uses a large number of Laplacian steps, n_steps.

Another point that led us away from a scalar Laplacians for velocity components is the question of what should be conserved, see the discussion in #26 and https://github.com/ocean-eddy-cpt/gcm-filters/issues/6#issuecomment-783625490. As @iangrooms points out above, it is natural to conserve a scalar quantity like salt via an integral over the full sphere. It is not natural to conserve the integral of a scalar component of velocity on the sphere. In other words, the global integral of a scalar component of velocity on the sphere is not a physically meaningful quantity. Instead, we want to conserve angular momentum. And @StephenGriffies's (and MOM5's) lateral friction operator satisfies this property.

I started to implement the vector Laplacian on a C-grid because I need this operator for an offline analysis of NW2 and MEKE anyway. If you, @arthurBarthe and @LaureZanna want to tackle the B-grid case, that would be awesome!

[StephenGriffies]

I originally did not think my input would be of use to this discussion on filter details. But given the strong overlap with diffusion and friction operators implemented in ocean models, it seems I do have stuff to add to the mix. Although I am not coding, I am happy to offer further input on this work either by looking at code and/or further discussions of the theory. So please keep me entrained if you have other related threads.

[Hallberg-NOAA]

The C-grid vector filter that is being discussed here would seem to me to be entirely analogous to what we are already doing in MOM6 in the operator that we use for the Laplacian viscosity, (For the code see https://github.com/NOAA-GFDL/MOM6/blob/dev/gfdl/src/parameterizations/lateral/MOM_hor_visc.F90 .) Although this file is large (2800 lines), most of that is about how the viscosity itself is specified or about the inclusion of biharmonic viscosities. This MOM6 viscous acceleration is implemented as the divergence of a stress tensor, and it includes all of the metric terms for use on a locally orthogonal grid (like a tripolar grid). Because of the form that is chosen, it is guaranteed to conserve a discrete form of angular momentum, and it includes a mask to prevent values from land from influencing the ocean.

If this MOM6 viscosity code turns out to be the template for what we want to do for the filter for vectors, we might want to treat the thicknesses as uniform, if it is the velocities or stresses and not the layer-integrated momentum that is being filtered. For a globally uniform length scale, the analog would be the viscous acceleration with a constant viscosity. The metric terms used here are derived from the basic grid metrics stored in the MOM6 ocean grid type defined in MOM_grid.F90 (https://github.com/NOAA-GFDL/MOM6/blob/dev/gfdl/src/core/MOM_grid.F90).

[NoraLoose]

Thanks @Hallberg-NOAA for pointing me to the relevant places in the MOM6 code! Mimicking the MOM6 implementation to define a vector Laplacian on a C-grid has been a very useful strategy.

I am currently writing some tests that check whether my vector Laplacian conserves angular momentum, as it should. I'm wondering if angular momentum conservation is part of the MOM6 testing framework, too?

[NoraLoose]

While working through the MOM6 MOM_hor_visc code, I noticed an inconsistency in how the biharmonic friction operator is implemented compared to what Griffies & Hallberg (2000) suggest:

• MOM6: div. ( Ah grad ( div . grad U ) );
• Griffies & Hallberg (2000): div. ( sqrt(Ah) grad ( div . (sqrt(Ah) grad U ) ).

The two options are equal if Ah is spatially constant. As Griffies & Hallberg (2000) note, the first operator is not guaranteed to dissipate kinetic energy (unless Ah is spatially constant), but the second operator is.

For the gcm-filters scalar (diffusion-based) biharmonic operator, we have been following the second recipe, and this is how I was going to implement the vector Laplacian too. But I am curious to hear if the MOM6 implementation has its own advantages. Any thoughts on this @adcroft, @StephenGriffies, @Hallberg-NOAA?

[StephenGriffies]

We recently had this question from Andrew Stewart for his implementation of biharmonic mixing in another model. I will defer to @Hallberg-NOAA for reasons why he did not follow the sqrt approach discussed in Griffies and Hallberg. I do not recall it being a conscious decision, perhaps only an oversight?

For the gcm-filters, I recommend following the Griffies and Hallberg approach since that indeed guarantees KE is dissipated.

[iangrooms]

We've got one vector Laplacian, and developing another in #106, so I'm closing this.