Optimized connectivity: Simultaneous face_edges and edge_nodes

[cmdupuis3]

Would close #1138, #1196

Related to #1180

Supercedes #1195

Overview

This set of changes optimizes face_edge, edge_node, and face_face connectivity. face_edge and edge_node connectivity are combined into one routine, while face_face is optimized stand-alone.

PR Checklist

General

• An issue is linked created and linked
• Add appropriate labels
• Filled out Overview and Expected Usage (if applicable) sections

Testing

• Adequate tests are created if there is new functionality
• Tests cover all possible logical paths in your function
• Tests are not too basic (such as simply calling a function and nothing else)

Documentation

• Docstrings have been added to all new functions
• Docstrings have updated with any function changes
• Internal functions have a preceding underscore (_) and have been added to docs/internal_api/index.rst

[cmdupuis3]

I think all the hard parts of the merge are done.

At this point, the only failures seems to be sorting issues.

[hongyuchen1030]

I think all the hard parts of the merge are done.

At this point, the only failures seems to be sorting issues.

@cmdupuis3 Thanks for your work, I will review it as soon as possible.

And do you have any idea why the CIs are all failing?

[cmdupuis3]

The CI is failing because the new algorithm returns the results in a different order.

I was thinking we can add some sorting mechanism, at least for legacy behavior. Phillip was evidently aware of this issue too. It depends on if you need to support that... Personally I'd be okay with just changing it, but I think you would be a better judge of the situation.

[cmdupuis3]

Altering njit and .values here runs the risk of conflicting with #1583. I'm thinking we can scope it so that .values cleanup in connectivity.py is allowed on this PR and everywhere except connectivity.py belongs to #1583.

[Sevans711]

Took a close look at this one just now. It does look like it is getting closer to being ready overall, I have a few small questions and I left some inline comments accordingly.

The main blockers for me would be a couple bigger questions which I do not think are resolved yet. These are: "why does this PR introduce lexicographic ordering?" and "do the changes here actually speed up n_nodes_per_face for large grids?" (See replies to older comment threads for more details).

[cmdupuis3]

pre-commit.ci autofix

[erogluorhan]

Please see below a few comments:

[erogluorhan]

Can you elaborate on this?

If the input grid has edges defined but no explicit edge-node connectivity, will we never be able to create that connectivity?

Even if that's the case, UXarray should still be able to construct connectivity that respect the edge ordering in the existing edges, and I believe this has been the case so far with our current connectivity construction.

[cmdupuis3]

So, what's happening here is that if a grid already comes with edge_node connectivity, that would be loaded. It would presumably be in some order, but we don't know what order that actually is. If _populate_edge_node_connectivity gets called on that Grid, it's going to reconstruct the edge_node connectivity in lexicographic order, which may not match the original order. So the idea is that if you have variables indexed by edge_node, recalculating edge_node connectivity could potentially scramble the index order.

I changed this message to hopefully be clearer about what's happening.

[erogluorhan]

That sounds good, but what is concerning to me is that this if-check will trigger even when there is only edge coords but not edge connectivities present in the grid, e.g. edge_lon, edge_lat. What do you think?

[cmdupuis3]

@erogluorhan Maybe if "edge_node_connectivity" in grid._ds would be more accurate?

[erogluorhan]

Yeah, I think that is safer

[rajeeja]

No, it's not reachable. All three callers check "edge_node_connectivity" not in grid._ds before calling, so a grid that already has edges from a file returns them straight from the property and never gets here.

[erogluorhan]

You're right, callers already handle that check.

@cmdupuis3 if you want to handle this check only here, that's okay, but remove the same check from all the other callers.

Also on a second thought, should we really raise a ValueError or just silently skip? I think the latter.

[Sevans711]

I don't think any other _populate... functions have internal logic like "if check: actually, don't populate anything". It would feel strange to me to add that here, it would not follow my expectations that "calling _populate...() should actually populate something". If there is some reason it can't actually run as expected, I would want it to crash (with ValueError if it is fundamentally impossible for the given inputs, or NotImplementedError if it is possible but would just need a cleverer algorithm which doesn't exist yet).

[erogluorhan]

We are talking about not populating again when there is already edge_node_connectivity present in the Grid (i.e. in order not to possibly break lexicographic order in this case). It can occur in several cases though, i.e. either when grid.edge_node_connectivity is accessed after the first time, or grid.face_edge_connectivity is being populated but edge_node_connectivity already exists.).

I don't agree with throwing an error, i.e. making code crash in that case doesn't make sense since the user has everything ready, and why should avoidance of re-populating an existing property crash the code, shouldn't it be okay with a warning instead?

If we don't want if-check within _populate...(), remove it since all of the callers are already doing this check and move the above message to where it is called (only for the case of simultaneous construction with face_edges I believe)

[cmdupuis3]

@rajeeja Good catch, after looking at the call sites, it looks unreachable to me. On digging deeper though, I found that where Philip left this branch, this block was included but commented out:

    # if (
    #     "edge_node_connectivity" not in grid._ds
    #     or "inverse_indices" not in grid._ds["edge_node_connectivity"].attrs
    # ):
    #     _populate_edge_node_connectivity(grid)

This would have been a case where you'd want to have the check on if "edge_node_connectivity" in grid._ds:. But part of the point of this branch is to remove the need for inverse_indices, so without that, the check becomes degenerate.

I think with that, it makes sense to just remove it.

@erogluorhan To me, it's about keeping the state of Grid consistent. Having edge_node_connectivity in a different state than the rest of the Grid object seems like it would get very confusing. But as Rajeev pointed out, this particular guard is unreachable anyway.

[Sevans711]

All of my prior comments have been resolved and no glaring issues popped out at me from a final read-through, so I'm happy to approve!

[rajeeja]

Going through this in pieces, will have a few more comments later today. Good PR overall and happy to be reviewing it.

[rajeeja]

@cmdupuis3 The new ValueError in _populate_edge_node_connectivity can't be reached from any caller, so it's untested and shows as uncovered. This MPAS mesh already loads with edge_node_connectivity in _ds, so it pins the guard without a new fixture.

Suggested change

	assert np.all(face_edges[face_idx, n_edges:] == INT_FILL_VALUE)
	assert np.all(face_edges[face_idx, n_edges:] == INT_FILL_VALUE)
	def test_connectivity_edge_node_refuses_to_overwrite(gridpath):
	"""Test that rebuilding edges on a grid that already has them is refused."""
	from uxarray.grid.connectivity import _populate_edge_node_connectivity
	# This mesh supplies verticesOnEdge, so the grid loads with edges already present
	uxgrid = ux.open_grid(gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc"))
	assert "edge_node_connectivity" in uxgrid._ds
	with pytest.raises(ValueError, match="already has"):
	_populate_edge_node_connectivity(uxgrid)

[rajeeja]

Take outCSne30.ug. Renumber its face_edge_connectivity edge ids, put on a fresh grid with
no edge_node_connectivity. Legal UGRID.

• Touch grid.n_edge. Main: 5400/5400 rows survive. This PR: 5400/5400 overwritten. No error,
  no warning.
• Cause: the PR builds both connectivities in one pass, so _populate_edge_node_connectivity
  now writes two arrays. The guard at connectivity.py:173 still checks one.
• Second write at connectivity.py:200 walks past it.
• The guard's reasoning is right. Constructed edges are lexicographic and won't match file
  order, which is what its own error message says. It just needs to cover both variables the
  function writes now.
• Main isn't correct here either, it keeps the file's array while the new
  edge_node_connectivity is lexicographic, so the two disagree silently.
• Untestable today: no mesh under test/meshfiles has face_edge_connectivity without
  edge_node_connectivity. Same hole as test_connectivity.py:176, one variable over.

Two smaller things in utils.py:529-536. The comment on MIN_ADAPTIVE_SORT_SIZE = 16 says a
bucket that size can't exceed the shift budget, which holds only because
MAX_SHIFTS_PER_EDGE = 8: max inversions are size(size-1)/2 against a budget of 8*size, so
size <= 17. That coupling isn't written down or asserted anywhere, and dropping MAX to 4
makes the comment false at size 9.

And test_connectivity_bucket_sort says the 500 bucket falls back to the heap sort, which it
does, but the test only checks the final order and insertion sort alone would produce the
same. Raise MAX_SHIFTS_PER_EDGE and the heap path goes uncovered with the test still green.
_heap_sort_bucket isn't imported by any test.

[cmdupuis3]

@rajeeja Wait, so if I understand correctly, this issue is more that the guards at call sites don't check both connectivities, so without having if "edge_node_connectivity" in grid._ds or "face_edge_connectivity" in grid._ds:, we still have the possibility of desynching between these two connectivity variables, right? It would just be a different mode of desynch than I was thinking of.

I'll see about consolidating the two sorting parameters into one and moving it to constants.py.