Conservative velocity regridding that conserves global transports

[taimoorsohail]

Hi!

I have been interested in computing the conserved velocity when moving from a curvilinear grid like a Tripolar grid to a rectilinear/lat-lon grid in Oceananigans.jl.

I think the correct approach would be to identify the intersection faces of the source grid that lie within the destination grid cell. We could then integrate u and v on these face to get a net u_transport_src and v_transport_src. We then rotate these two components to be orthogonal to the destination grid and compute the new u_transport_dst and v_transport_dst there. Finally, we divide by the destination grid cell face area to get u_dst and v_dst.

The key issue here is identifying the intersection faces to integrate over, particularly as most are partial grid cell faces. I think I have an idea that I can add to the Oceananigans extension.

cc @simone-silvestri @navidcy @glwagner

[asinghvi17]

The way I had imagined doing velocity regridding would be that you transform all velocities to be in unit spherical coordinates, then use the usual regridder, and then perform the transform to cell-local tangent space. I'm by no means an expert here though, so is there some way in which that's not sufficient? Happy to chat more.

[taimoorsohail]

So the key here is transport conservation for velocities which requires conserving the face-flux transports (i.e. v*dA_y; u*dA_x).

The approach you propose rotates vectors into unit-spheres, then regrids like a scalar, so it conserves cell area averages, not face transports.

In the PR, I integrated source face transport over overlaps with the destination grid (in both x- and y-), rotated that face transport onto destination faces, then divided by destination face area to get u/v. This approach exactly conserves the net eastward and northward transport. The u and v may change for a given grid cell, but the total transport stays exactly the same (within 1e-14m/s). This is essential for velocity-based diagnostics like stream functions which require global closure of transports.

I hope this makes sense!

[navidcy]

Just noting that in Oceananigans, because of the Arakawa C grid, the u and v velocities live on different points on the stencil.

https://clima.github.io/OceananigansDocumentation/stable/appendix/staggered_grid

[asinghvi17]

hmm, I think I see now. The way this is done in the PR isn't super efficient, and I'd like to avoid depending on LibGEOS. So I'd propose we sketch out the APIs here - what the regridder needs to look like, what we need from input data that is being regridded, etc. - and then I can propose a method that works?

This probably also wants some guarantee of an order of points in each cell for at least Oceananigans, so we know from the source grid what lines are u and what lines are v.

[navidcy]

This probably also wants some guarantee of an order of points in each cell for at least Oceananigans, so we know from the source grid what lines are u and what lines are v.

Not sure I understand what you mean here?

[navidcy]

I'm not sure what you mean also with "what the regridder needs to look like". But the main difference for regridding a vector field is that we need to regrid both u and v together.

So the regridder needs to take both e.g. u_src and v_src and u_dst and v_dst as inputs where u_src and u_dst are Oceananigans XFaceFields that live on the source and destination grids respectively (and v_src, v_dst the same but YFaceFields).

[navidcy]

After the regridder is constructed, perhaps I suggest that we add a method so that we regrid via

function ConservativeRegridding.regrid!(dst::Tuple{XFaceField, YFaceField}, Regridder::VelocityRegridder, src::Tuple{XFaceField, YFaceField})
    ...
end

And we apply the regridding like:

R = VelocityRegridder(....)

ConservativeRegridding.regrid!((dst_u, dst_v), R, (src_u, src_v))

[navidcy]

The above should work independently for every vertical level, so we could launch a kernel in k perhaps?

[taimoorsohail]

The most complex thing here is to compute the intersection faces (and associated weights) between the source and destination grids. I used LibGEOS to do this, hence the dependency. If we can come up with a way to compute line intersections without it I would be happy if we dropped the dependency.