Tight axis limit for planar pediodic meshes

[andrewdnolan]

Per some discussion in E3SM-Project/polaris#256, we should be explicitly setting the axis limits for planar periodic plots.

When repeating of periodic patches across boundaries is supported this will produce a actual periodic plot.

[andrewdnolan]

Testing

import mosaic
import matplotlib.pyplot as plt

def test_periodic(descriptor, ds, loc):
    fig, ax = plt.subplots()

    mosaic.polypcolor(ax, descriptor, ds.indexToCellID * 0.0, ec='k', cmap='binary')
    mosaic.polypcolor(ax, descriptor, ds[f"indexTo{loc}ID"], alpha=0.5, ec='tab:orange')

    ax.scatter(ds[f"x{loc}"], ds[f"y{loc}"], color="tab:blue")

    ax.set_title(f"{loc} Patches")
    
    fig.savefig(f"{loc}.png", dpi=300, bbox_inches='tight')

ds = mosaic.datasets.open_dataset("doubly_periodic_4x4")

descriptor = mosaic.Descriptor(ds)

for loc in ["Cell", "Edge", "Vertex"]:
    test_periodic(descriptor, ds, loc)

plt.show()

Cell patches:

Edge patches:

Vertex patches:

[xylar]

@andrewdnolan, this looks great! Would it make sense to also address the edge polygon that's on the wrong side of the periodic boundary here or do you want a separate PR for that?

[xylar]

I think the current strange xEdge locations come from the MPAS mesh creator. I don't see anywhere (maybe I missed it) that mosaic is trying to correct xEdge (internally xCell) for periodic meshes, only its internal xVertex, right? I think it's currently trusting the input mesh to have a reasonable xEdge (or xCell or xVertex) and it turns out it shouldn't.

[xylar]

From the results you're showing, I think this probably should instead be yVertex here instead of yEdge. And the same for xVertex and xEdge (which happen to have the same min). For now, that will cut off more of the top of cell, edge and vertex patches but eventually it will look right and it will start the figure at 0,0 as I think we would want. That way, the x and y periods are visible as the limits of the axes.

[andrewdnolan]

I agree, that it should be yVertex instead of yEdge. If you take a look at where mesh variables are defined below, we'll actually need to use yVertex.max() - y_period. That's fine for now, but I'm just a little curious how generalizable this is? Are the only periodic meshes ever generated by the make_planar_hex_mesh function? If that's the case, then I guess we can be confident this will generalize.

[andrewdnolan]

@xylar Following on my response you your comment above. Is there ever any chance of the meshes being rotated 90 degrees? (i.e. so that the left most edges are not parallel with the y-axis).

I guess the way to be most general would be to use the most extreme vertex position, (i.e. something like this):

y_min_edge = float(descriptor.ds.yEdge.min())
y_min_vertex = float(descriptor.ds.yVertex.min())

if y_min_vertex > y_min_edge: 
    ymax = float(descriptor.ds.yVertex.max())
    ymin = ymax - descriptor.y_period 
else: 
    ymin = float(descriptor.ds.yVertex.min())
    ymax = ymin + descriptor.y_period 

using the xVertex something like whats above would also probably catch if the domain was ever rotated 90 degrees.

[andrewdnolan]

@xylar Re the xEdge location outside of the proper periodic, I'd vote to punt on that in this PR.

I did a little minimal testing and I think that misplaced edge can be implicitly corrected when we mirror patches across the periodic boundaries. There's some minor details to work out, but I think I have a pretty good idea of how we can implement this for planar periodic meshes. If this misplaced edge can get fixed with the mirroring, I'd consider that preferable than having an additional method/function to deal with this misplaced edge (given that mosaic doesn't really consider it misplaced).

[xylar]

I guess the way to be most general would be to use the most extreme vertex position

I would think taking the minimum vertex position would always be safe but I could be wrong. Your approach likely won't cause trouble.

[xylar]

We don't rotate periodic meshes in practice but it could be done so, yes, you might as well handle it.

[xylar]

I'm definitely fine with you fixing the edges on the wrong side of the periodic boundary "properly" by mirroring if you think that's not too arduous.

[andrewdnolan]

@xylar The reason I'm opting for the more complicated approach, in the most recent commit (87fbd38), is because that still allows us to keep all patches agnostic of each other.

While I agree with:

I would think taking the minimum vertex position would always be safe but I could be wrong.

If you take a look at the attached plot, you'll notice the scattered vertex positions in purple don't appear at the the minimum vertex position of the corrected patches. If we just used the minimum vertex position from the descriptor.ds.yVertex array, we end up getting the incorrect axis limits.

We could get the minimum from the corrected cell patches, but that means we'd have to create the cell patches in order to set the axis limits of the edge and/or vertex patches, which goes against the lazy loading approach we've implemented up until this point. I'm pretty confident that what's implemented in the most recent commit will correctly set the axis limits, while keeping patches independent of one another.

[xylar]

Yes, I agree, what you have implemented is better and makes sense to me.