Fix planarIntersect formula in MPAS mesh converter

[mwarusz]

The new formula is from: https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line

[xylar]

@mwarusz, have you checked whether a similar fix is needed for spherical meshes?

[xylar]

A question to @mwarusz and reviewers: is it not the case the correct edge location halfway between vertices (see #565 (comment)) rather than at this intersection?

Perhaps the intersection formula might be correct but even if so, is it actually the case that we want the edge to be located halfway between cell centers along a great circle rather than at the intersection?

[xylar]

Finally, currently the edge location for edges with only one neighboring cell are halfway between vertices:

MPAS-Tools/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp

Line 1278 in d29c4f2

edge_loc = (vert_loc1 + vert_loc2) * 0.5;

I don't see a good alternative to this but wanted to bring it up.

[favba]

A question to @mwarusz and reviewers: is it not the case the correct edge location halfway between vertices (see #565 (comment)) rather than at this intersection?

Perhaps the intersection formula might be correct but even if so, is it actually the case that we want the edge to be located halfway between cell centers along a great circle rather than at the intersection?

@xylar , The mesh specification is indeed a bit confusing regarding that.
In the second paragraph of chapter 2 it says:

As shown in Figure 2.1, a line segment that connects two primal mesh cell centers is uniquely
associated with a line segment that connects two dual mesh cell centers. We assume that these
two line segments cross and the point of intersection is labeled as x_e .

Which is what me and @mwarusz have been suggesting.

But later, in the first paragraph of Chapter 3 it also says:

In MPAS, cells
are nominally located at the Voronoi generating points, which, for centroidal Voronoi tessellations,
are the mass centoids of the Voronoi cells with respect to a density function, and edges are nominally
located at the midpoints of edges.

Which suggests your point of view. But I take that that "nominally" in the text means it isn't truly there...

I'll also note that my advisor published a paper in 2016 where, one of the things analyzed, is precisely the current edge positioning of MPAS (at intersections) vs a modified method that uses the edge midpoint.

[mwarusz]

@mwarusz, have you checked whether a similar fix is needed for spherical meshes?

I haven't checked the formula, but it seems to be doing the right thing based on numerical testing. Note that the current behavior of the spherical version and the planar version on distorted meshes is opposite with respect to edge point placement.

A question to @mwarusz and reviewers: is it not the case the correct edge location halfway between vertices (see #565 (comment)) rather than at this intersection?

Perhaps the intersection formula might be correct but even if so, is it actually the case that we want the edge to be located halfway between cell centers along a great circle rather than at the intersection?

I believe this recent paper https://www.sciencedirect.com/science/article/abs/pii/S146350032100161X
from your group discusses the difference between the two approaches and concludes that what this PR does is better. Let's wait to hear what @mark-petersen has to say.

[xylar]

@matthewhoffman and @trhille, MALI is likely to be the component most affected by this change, since you use irregular, planar meshes more often than the other components do.

[xylar]

@mgduda, this is one you should definitely weigh in on if you have time.

[xylar]

I tested this on our Compass ocean pr test suite. It isn't bit-for-bit with the version using the previous, incorrect intersection calculation but differences were at machine precision, which is what I would have expected. All of our planar meshes are regular so the change in edge position would be expected to be very small.

[xylar]

I will make an mpas_tools release candidate with this fix in it to help others test it.

[xylar]

A package with this change can be installed with:

conda create -y -n test -c conda-forge/label/mpas_tools_dev mpas_tools=0.34.0rc1

or in an existing environment with:

conda install -y -c conda-forge/label/mpas_tools_dev mpas_tools=0.34.0rc1

[mark-petersen]

Thanks. I agree with these changes based on @xylar's explanation and results.

[xylar]

More testing

I ran @favba's mesh through the converter with this update.

Before the fix, the edges and cells around the 13th vertex looked like this:

With this fix, it looks as expected:

[xylar]

Approved based on my compass testing and plotting.

[mgduda]

For MPAS-Atmsophere, we use the TRiSK placement of edge locations (i.e., at the midpoint of triangle faces), consistent with the changes in this PR.

[proteanplanet]

This is to accept the change on the planar mesh, and to confirm that for spherical meshes, the edge point falls at the intersection of the great circles between vertices and cell centers. I have confirmed this graphically for the MPAS-Ocean/SeaIce EC30to60E2r2 mesh, centered here for the Pacific 14.5˚N 166.4˚E on a 0.75˚ horizon:

Mesh lines on this plot are great circles at 1 NM resolution. Edge points as defined in mesh files are plotted as blue dots.

[xylar]

@proteanplanet, thank you so much, that's reassuring!

[matthewhoffman]

Based on the discussion and inspection of the changes, this is convincing. @trhille and I will test this with a few of our mesh generation workflows before approving to ensure it doesn't generate any surprises for us.

[matthewhoffman]

I ran a few of our compass mesh generation test cases for nonuniform meshes generated via JIGSAW using the rc mpas-tools package that @xylar conveniently set up.

As expected, all meshes generated without incident and look identical in Paraview, even with the vtk converter. At first I expected to see very slightly different edge locations, but I believe the vtk polygons don't actually use xEdge, yEdge, and those are the only fields modified by this PR. I did confirm that xEdge, yEdge differ slightly before and after this PR, and that xCell, yCell, xVertex, yVertex, dcEdge, and dvEdge all are identical before and after. Given that we don't use xEdge,yEdge anywhere in MALI (other than the mpas_li_velocity_simple.F module which hasn't been used in years), I don't anticipate this PR to actually impact simulations in any way.

So to summarize, no impact on MALI simulation, no impact on common visualization tools, and clearly a correct change - I'm happy to approve and let this move forward.

[xylar]

@matthewhoffman, that's wonderful! I expected MALI to be impacted. If not, I think we can merge this without concerns.

@mwarusz, I really appreciate your diligent and thoughtful work on this.

@favba, I deeply appreciate you bringing this up and engaging with us long enough for us to both appreciate that it was a real, significant problem and to find a solution we're happy with.

[xylar]

I will make a new MPAS-Tools release soon with this fix. I just need to check if anything else should be included.