mesh: compute 1D inter-process facets from the vertex IndexMap

[jhale]

Topology's tdim==1 constructor hack tried to derive inter-process
facets (vertices, for a 1D mesh) via mesh::compute_entities(dim=0),
but that function has always returned immediately for dim==0 since
vertices are set up directly rather than computed on demand. The
pushed (always-empty) vector satisfied interprocess_facets()'s "not
computed" guard, so it silently returned an empty list for every 1D
mesh instead of the true inter-process vertices, misclassifying
partition-boundary vertices as exterior facets under GhostMode::none.

Replace the hack with compute_interprocess_vertices, which reads the
answer directly from the vertex IndexMap's own ownership/ghosting
data instead of recomputing entity sharing: every vertex needed by a
locally owned cell is already resolved there (owned or ghost),
independently of GhostMode, so no extra communication round is
needed. Drops the now-unused num_threads parameter from Topology's
constructor.

Co-Authored-By: Claude Sonnet 5 noreply@anthropic.com

---

_{Stack created with GitHub Stacks CLI • Give Feedback 💬}

[jorgensd]

@jhale can you check that this is correct for t joints or even larger joints? Ie sometimes where a facet is shared between 3 to N cells?

[jorgensd]

The compute_interprocess_vertices discussion made me look at index_to_dest_ranks, where the observation that the IndexMap already carries the communication pattern is not currently exploited.

Therefore I asked CLAUDE (OPUS 5) to see if we could exploit this rewrite in index_to_dest_ranks. The following is a summary from my discussion with CLAUDE:

Follow-up idea (not blocking): IndexMap::index_to_dest_ranks could be rewritten along the same lines**

I instrumented a copy of the current implementation (verified to return identical data/offsets) and timed it on the vertex map of a 32³ tetrahedral cube, GhostMode::shared_facet, 8 ranks, max over ranks, 20 reps:

index_to_dest_ranks (full)                   : 0.211 ms
  consensus NBX (dest discovery)             : 0.076 ms
  part 1: ghosts -> owner exchange           : 0.099 ms
  part 2: build send buffer                  : 0.013 ms
  part 2: comm (create comm + alltoallv)     : 0.053 ms
  part 2: unpack (2 sorts + binary searches) : 0.068 ms

Scatterer construction (same exchange as part 1)   : 0.055 ms
owner-side rows rebuilt from Scatterer data (local) : 0.003 ms

Three changes, in order of payoff:

1. Drop the consensus round. The function calls compute_graph_edges_nbx to discover dest, but the constructor already computed exactly that into _dest via build_src_dest. I checked on 4 and 8 ranks that the rediscovered list equals dest() on every rank. A caller-supplied dest is only verified under #ifndef NDEBUG, but a superset merely contributes empty blocks to the neighbourhood exchange, and a subset is a map that is already broken for Scatterer, shared_indices and every scatter operation, since they all trust _dest. So this is pure subtraction, and it removes a global MPI_Ibarrier-based synchronisation point that will cost more than 0.076 ms at high rank counts. It also makes the tag parameter dead, so the parameter and its @note in IndexMap.h can go (public API change).

2. Part 1 is the Scatterer constructor. Both send ghost global indices to their owners over the same src/dest neighbourhood, and Scatterer::local_indices_block() together with _sizes_local is exactly the (owned local index, ghosting rank) list that part 1 rebuilds with an int64 pair sort. Given that data the owner-side rows are a counting sort: 0.003 ms against 0.099 ms, and one MPI_Dist_graph_create_adjacent disappears. This needs Scatterer to expose its per-neighbour sizes/displacements, or the regrouping to move inside Scatterer/IndexMap.

3. Part 2: send in forward-scatter order. Today each owned index i with n_i ghosting ranks sends (global index, rank) int64 pairs to each of them, 2·Σn_i² int64 in total, and the receiver sorts its ghosts, binary-searches every received index, then sorts again. In scatter order the buffer position already identifies the ghost through remote_indices_block(), so the global index need not travel at all: int32 counts plus int32 rank lists, and the unpack becomes one linear pass. Payload on the mesh above: 146 KiB → 67 KiB, approaching 4× as n_i grows (corner vertices, high rank counts).

Extrapolating from the section timings, and without having prototyped it, that is roughly 0.21 → 0.07–0.09 ms, so 2–3×.

Independently of the above, get_local_indexing calls index_to_dest_ranks once per entity dimension on the same vertex map, so a 3D mesh that creates both edges and facets computes the identical result twice (more with mixed topology). Hoisting the call, or caching the result on the map, halves the cost before anything else changes. The other callers are refinement/interval.h and refinement/plaza.h.

Caveats on the numbers: shared-memory MPI on a 10-core machine at ≤ 8 ranks, so the balance between local work and communication will look different at scale; the per-section maxima land on different ranks and therefore do not sum to the total.

[jorgensd]

This is a very cool way of doing ghosting, which I believe could be used to extend the number of layers one would search outwards in the dual graph! @IgorBaratta

[schnellerhase]

Which use cases does this fix? For the here tested branching manifolds, there is no general concept of interior or exterior facet.

[jhale]

Which use cases does this fix? For the here tested branching manifolds, there is no general concept of interior or exterior facet.

The PR fixes the computation interior and exterior facets with GhostMode.None for one-dimensional non-branching meshes (see e.g. the super simple included test with just a UnitInterval, nothing fancy). Don't get overly distracted by non-manifold meshes, this is just test harnessing to check the definition (more below!) is actually applied in more complex cases.

Before the fix, interprocess_facets() always returned empty, because the code path it relied on compute_entities(dim=0) which returns empty. With GhostMode::none, a vertex at a partition boundary has only one cell attached on its own rank, so it "looks" exterior to itself unless interprocess_facets() says something otherwise. This all mysteriously worked with GhostMode::SharedFacet because of the automatic duplication of ghost geometry data on interior facets at partition boundaries.

The definition of interior and exterior is:

Let $f$ be a facet (a $(d-1)$-dimensional entity of a mesh of topological dimension $d$), and let $c(f)$ denote the number of cells owned globally that have $f$ as a boundary entity.

Exterior facet: $c(f) = 1$.
Interior facet: $c(f) \geq 2$.

This definition already exists, from my understanding of the existing code.

[schnellerhase]

Thanks, should this not be managed in the create_entities special case handling, rather than in the constructor?

1. Computing interprocess facets data in the constructor forces computational cost which is not always required,
2. We expect any code path to rely on facet information to call a mesh.topology.create_entities(tdim-1) at some point anyhow, so this should be consistent for tdim=1 (for tdim=1 we would have the special case that vertices are always available but interprocess facet data might not be available yet)

One other point, that I need to think more on: if the shared facet case worked before without this additional work path than there must be a no-op short cut to this computation at least in this case, which we should retain.

[jhale]

Thanks, should this not be managed in the create_entities special case handling, rather than in the constructor?

Yes.

1. Computing interprocess facets data in the constructor forces computational cost which is not always required,

Indeed, an easy fix, I will implement it through lazy computation.

2. We expect any code path to rely on facet information to call a mesh.topology.create_entities(tdim-1) at some point anyhow, so this should be consistent for tdim=1 (for tdim=1 we would have the special case that vertices are always available but interprocess facet data might not be available yet)

Agreed.

One other point, that I need to think more on: if the shared facet case worked before without this additional work path than there must be a no-op short cut to this computation at least in this case, which we should retain.

It can be gated on the same max_cell_facet_links parameter that you introduced, and I guess GhostMode too which is lucky because @garth-wells fixed that very recently! But as it requires a bit of routing through, so I will do it on a follow up branch.

Edit: The old GhostMode::shared_facet route was a bit mysterious, but perfectly legitimate, and this is just a fix for GhostMode::none.

[garth-wells]

Comments inline.

No chance that we can get tdim==1 to follow the same code path as tdim>1?

[garth-wells]

It's not clear what an 'inter-process' vertex is.

[jhale]

It's described in the docstring:

/// Matches the definition used by mesh::compute_entities for dim >= 1:
/// an entity is inter-process if cells owned by two or more ranks are
/// attached to it.

And is consistent with the interprocess_facets which is in the public API.

[garth-wells]

Use transforms in cases like this is harder to read than a plan loop.

[garth-wells]

Add brace for readability.