Consistency in defining Laplacian operators

[ElizabethYankovsky]

Hello everyone,

I and @arthurBarthe have been working on the kernel for MOM5, and after speaking with @NoraLoose want to make sure that the Laplacian operators we implement are consistent across the various kernels.

In MOM5 (B-grid, uses a northeast convention), the grid variables that we use to construct the Laplacian are:

• dxt, dyt: x- and y-spacings of the tracer cells
• dxu, dyu: x- and y-spacings of the velocity cells
• dau, dat: surface areas of velocity and tracer cells (respectively).
  Please see Figure 9.7 for details on the MOM5 grid: https://mom-ocean.github.io/assets/pdfs/MOM5_manual.pdf

Our Laplacian operator on a field a defined at velocity points is constructed by first computing first derivatives fx and fy:
fx(i,j) = (a(i+1,j)-a(i,j))*1/(dxt(i+1,j)); fy(i,j) = (a(i,j+1)-a(i,j))*1/(dyt(i,j+1))
We then compute LAP_U (Laplacian at u points) as:
LAP_U(i,j) = (dyu(i,j)*fx(i,j) - dyu(i-1,j)*fx(i-1,j))*1/(dau(i,j) + (dxu(i,j)*fy(i,j) - dxu(i,j-1)*fy(i,j-1))*1/(dau(i,j)

Is this in line with what others are implementing? We based our definition on the LAP_T operator defined in the MOM5 https://github.com/mom-ocean/MOM5/blob/master/src/mom5/ocean_core/ocean_operators.F90, but are open to discussing other forms. In particular, @NoraLoose is using a flux approach in order to maintain conservation properties.

[ElizabethYankovsky]

Also, for reference here is a schematic of the NE convention B-grid on which the above is based:

[NoraLoose]

@ElizabethYankovsky, I think you are essentially using the same flux approach as me in #23 - the first half of your Laplacian would be the "uflux" (sum of fluxes across western and eastern cell edge), and the second half of your Laplacian would be the "vflux" (sum of fluxes across northern and southern cell edge).

With the following minor modifications, your code would satisfy conservation properties, too (and could easily handle a no-flux boundary condition at the continental boundaries):

fx(i,j) = (a(i+1,j)-a(i,j))*2/(dxt(i+1,j)+dxt(i+1,j+1)) 
fy(i,j) = (a(i,j+1)-a(i,j))*2/(dyt(i,j+1)+dyt(i+1,j+1))

fx(i,j) and fy(i,j) would be the fluxes across the eastern and northern cell edge of the U(i,j) cell, respectively.

LAP_U(i,j) = - ( (dyu(i,j)+dyu(i+1,j))/2*fx(i,j) - (dyu(i-1,j)+dyu(i,j))/2*fx(i-1,j) ) * 1/dau(i,j) 
             - ( (dxu(i,j)+dxu(i,j+1))/2*fy(i,j) - (dxu(i,j-1)+dxu(i,j))/2*fy(i,j-1) ) * 1/dau(i,j)

[NoraLoose]

I agree with @ElizabethYankovsky: I think it should be clearly stated (somewhere obvious in this repo) which properties the Laplacian implementations should satisfy. As per my understanding, our design method crucially uses the property

Integral(Laplacian) = 0,

if you take the integral over the full oceanic domain. This property is guaranteed for finite volume implementations of the Laplacian with no-flux boundary condition at the continental boundaries. @iangrooms, please correct me if I'm wrong.

[ElizabethYankovsky]

Thanks @NoraLoose! I'll modify the Laplacian to be in the form you showed above. When you perform the filtering on the vorticity (rather than velocity) do you use a separate Laplacian specifically for t-cell fields? It appears that vorticity is defined at t-cell centers rather than u-cell centers in a B-grid model.

[NoraLoose]

When you perform the filtering on the vorticity (rather than velocity) do you use a separate Laplacian specifically for t-cell fields? It appears that vorticity is defined at t-cell centers rather than u-cell centers in a B-grid model.

Yes, filtering a field defined on t-cells (rather than u-cells) requires a different Laplacian: one with updated (staggered) grid information. However, the underlying Laplacian operation is always the same, equal to the generic Laplacian IrregularCartesianLaplacianWithLandMask in #23 :

class IrregularCartesianLaplacianWithLandMask(BaseLaplacian):
    """̵Laplacian for irregularly spaced Cartesian grids with land mask.

    Attributes
    ----------
    wet_mask: Mask array, 1 for ocean, 0 for land
    dxw: x-spacing centered at western cell edge
    dyw: y-spacing centered at western cell edge
    dxs: x-spacing centered at southern cell edge
    dys: y-spacing centered at southern cell edge
    area: cell area
    """


    wet_mask: ArrayType
    dxw: ArrayType
    dyw: ArrayType
    dxs: ArrayType
    dys: ArrayType
    area: ArrayType


    def __post_init__(self):
        np = get_array_module(self.wet_mask)


        self.w_wet_mask = self.wet_mask * np.roll(self.wet_mask, -1, axis=-1)
        self.s_wet_mask = self.wet_mask * np.roll(self.wet_mask, -1, axis=-2)


    def __call__(self, field: ArrayType):
        np = get_array_module(field)


        out = np.nan_to_num(field)


        wflux = (
            (out - np.roll(out, -1, axis=-1)) / self.dxw * self.dyw
        )  # flux across western cell edge
        sflux = (
            (out - np.roll(out, -1, axis=-2)) / self.dys * self.dxs
        )  # flux across southern cell edge


        wflux = wflux * self.w_wet_mask  # no-flux boundary condition
        sflux = sflux * self.s_wet_mask  # no-flux boundary condition


        out = np.roll(wflux, 1, axis=-1) - wflux + np.roll(sflux, 1, axis=-2) - sflux


        out = out / self.area
        return out

This Laplacian applies, no matter if one uses a B- or C-grid, or if one filters a field defined on t-cells or u-cells or v-cells. The only difference lies in the specification of the grid information, i.e., wet_mask, dxw, dyw, dxs, dys, area. For instance, your LAP_U code above is identical to the generic Laplacian IrregularCartesianLaplacianWithLandMask if one sets area=dau,

dxw =  0.5 * (np.roll(dxt, 1, axis=1) + np.roll(np.roll(dxt, 1, axis=1), 1, axis=0)),

and similarly for dyw, dxs, dys (assuming that your coordinates have the order (j,i)).

This is why in #23 (comment), I started advocating for using one generic Laplacian (e.g., IrregularCartesianLaplacianWithLandMask) as a base Laplacian, which the model-specific Laplacians could inherit from. If dxw, dyw, dxs, dys are not directly available from model output (as in your case), their calculation could happen in the post_init method for each model-specific Laplacian. This would have the following advantages:

• We would avoid writing redundant code by taking advantage of class inheritance.
• We would have a faster implementation: a Laplacian that uses a pre-calculated dxw will be faster than a Laplacian that computes 0.5 * (np.roll(dxt, 1, axis=1) + np.roll(np.roll(dxt, 1, axis=1), 1, axis=0)) in every iteration.

Of course this is just an idea, and I am very happy to hear and discuss alternative suggestions.

[ElizabethYankovsky]

@NoraLoose, I'm rewriting the Laplacian as you have here and just wanted to check with you on a detail. I think that the negative signs shouldn't be there:
LAP_U(i,j) = - ( (dyu(i,j)+dyu(i+1,j))/2*fx(i,j) - (dyu(i-1,j)+dyu(i,j))/2*fx(i-1,j) ) * 1/dau(i,j) - ( (dxu(i,j)+dxu(i,j+1))/2*fy(i,j) - (dxu(i,j-1)+dxu(i,j))/2*fy(i,j-1) ) * 1/dau(i,j)
Should instead be:
LAP_U(i,j) = ( (dyu(i,j)+dyu(i+1,j))/2*fx(i,j) - (dyu(i-1,j)+dyu(i,j))/2*fx(i-1,j) ) * 1/dau(i,j) + ( (dxu(i,j)+dxu(i,j+1))/2*fy(i,j) - (dxu(i,j-1)+dxu(i,j))/2*fy(i,j-1) ) * 1/dau(i,j)

This is coming from the fact that in the basic version we have: (dyu[i,j]*fx[i,j]-dyu[i-1,j]*fx[i-1,j])
Do you agree?

[NoraLoose]

I wrote the finite volume discretization of the Laplacian here:

...

The 4 terms in your LAP_U match exactly the 4 terms on the right hand side of (1) (modulo devision by area), in the same order. So yes, your corrected version (without the negative signs) is correct. 👍