Optimized connectivity: Simultaneous face_edges and edge_nodes - #1560
Optimized connectivity: Simultaneous face_edges and edge_nodes#1560cmdupuis3 wants to merge 19 commits into
Conversation
|
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? |
|
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. |
| pass | ||
|
|
||
| edge_nodes, face_edges = _build_edge_node_connectivity( | ||
| grid.face_node_connectivity.values, grid.n_nodes_per_face.values |
There was a problem hiding this comment.
Please think about if you can get rid of .values calls to be chunked-array compatible.
There was a problem hiding this comment.
On it, I think Philip was trying to scalarize here but it introduced this regression.
There was a problem hiding this comment.
So the best I could come up with is
face_nodes, n_nodes_per_face = dask.compute(
grid.face_node_connectivity.data, grid.n_nodes_per_face.data
)
edge_nodes, face_edges = _build_edge_node_connectivity(
face_nodes, n_nodes_per_face, grid.n_node
)...but I'm unclear on if we're committing to dask or not. Other parts of the repo imply that dask is still sort of an optional load, and doing this would pretty much require a dask import for most grids, unless we have some kind of numpy/dask conditional.
There was a problem hiding this comment.
Alright, I hacked it a bit to get this:
computed = xr.Dataset(
{
"face_nodes": grid.face_node_connectivity.variable,
"n_nodes_per_face": grid.n_nodes_per_face.variable,
}
).compute()
edge_nodes, face_edges = _build_edge_node_connectivity(
computed.face_nodes.data, computed.n_nodes_per_face.data, grid.n_node
)There was a problem hiding this comment.
Question for @erogluorhan, does uxarray have any plans to support chunked grids? I was under the impression that data variables can be chunked, but that it is always safe to assume the grid itself can be fully loaded into memory. Is that incorrect?
That said, using more dask-compatible syntax like this doesn't seem like it would have a downside (no need to block merging on this question).
The optimized edge builder deduped half edges with a numba hash map, which
numbered edges in first-encounter order. Edges had previously been numbered
lexicographically by their (min_node, max_node) pair, as a side effect of the
np.unique(..., axis=0) the hash map replaced.
Global edge index is a public identity: it indexes edge_lon/edge_lat, edge
centered data variables, and edge_node_distances, so renumbering silently
re-pairs user data with different physical edges. It also broke the five
TestQuadHexagon connectivity tests, which assert on edge_node, face_edge,
node_edge, edge_face and face_face -- all the same renumbering cascading
through the derived connectivities.
Sort as the dedup mechanism instead of hashing. Node indices are dense
integers in [0, n_node), so a counting sort buckets the half edges by their
first node without any comparisons, and sorting each bucket by its second node
leaves the duplicates adjacent -- the dedup then falls out of the same walk.
Buckets hold one entry per edge incident to a node, so on a real mesh they are
tiny (node degree, typically under ten) and an insertion sort finishes them.
A bucket above MAX_INSERTION_SORT_SIZE is heap sorted so that a degenerate
mesh cannot degrade the build quadratically; np.argsort is deliberately not
used there, as numba's implementation degrades badly on structured input.
Half edges are identified throughout by their flat face_node_connectivity
index, which is also the face_edge_connectivity slot they are written back to,
so the sort needs a single permutation array and no mapping back.
This is faster and leaner than the hash map it replaces. On a synthetic one
million face quad mesh, measured by peak RSS rather than tracemalloc, which
does not observe numba's typed dict allocations:
dict build 397.3 ms 239.5 MB
bucket sort 85.1 ms 91.6 MB
The five legacy tests now pass unchanged. Adds order invariant coverage for
the canonical ordering and the face_edge positional contract, plus high degree
nodes either side of the insertion sort threshold.
Also casts n_nodes_per_face back to INT_DTYPE, so that the builder is not
compiled a second time for int64, and restores a blank line dropped between
two top level functions.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have improved:
Benchmarks that have stayed the same:
|
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
Testing
Documentation
_) and have been added todocs/internal_api/index.rst