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 https://github.com/MPAS-Dev/MPAS-Tools/issues/565#issuecomment-2082590736) 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: https://github.com/MPAS-Dev/MPAS-Tools/blob/d29c4f218dbd9a9ec4b3bc5de8e797fb1c0e1d8a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp#L1278 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

[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]

@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.

[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.