Skip to content

Fix quadratic cost in adaptive mesh refinement's coarse/fine cell-map construction - #5326

Merged
pbrubeck merged 4 commits into
mainfrom
pbrubeck/adaptive-optimise
Aug 6, 2026
Merged

Fix quadratic cost in adaptive mesh refinement's coarse/fine cell-map construction#5326
pbrubeck merged 4 commits into
mainfrom
pbrubeck/adaptive-optimise

Conversation

@pbrubeck

@pbrubeck pbrubeck commented Aug 5, 2026

Copy link
Copy Markdown
Contributor

Summary

adaptive_parent_child_cell_maps resolved each fine cell's coarse parent with PETSc's DMLabelGetValue, called once per fine cell. DMLabelGetValue should resolve a point in O(1); instead it scans every stratum of the label linearly until it finds the one containing the queried point. set_adaptive_parent_label gives every coarse cell its own stratum value, so on an adaptively refined mesh the per-fine-cell loop was effectively O(nfine · ncoarse), not O(nfine). Profiling adaptive_multigrid.py under -log_view showed refine_marked_elements dominated not by PETSc's refine_sbr mesh transform itself, but by time attributed to this uninstrumented loop.

adaptive_parent_child_cell_maps now walks the label by stratum instead, with DMLabelGetStratumSize/DMLabelGetStratumIS, one call per coarse cell. That uses PETSc's O(1) value→stratum hash map and touches each fine cell exactly once, so the loop is O(nfine + ncoarse). Output is unchanged — cross-checked bit-for-bit against the original point-wise result on refined netgen meshes.

On the AdaptiveMeshHierarchy demo (first 10 Dörfler-marked refinement levels of an L-shaped netgen mesh), adaptive_parent_child_cell_maps's share of time inside refine_marked_elements drops from 59% to 8%, and refine_marked_elements's own total drops by 55%. The dominant cost left inside refine_marked_elements is now legitimately PETSc's own mesh-refinement transform.

Collateral changes

_adapt_marked_cells and refine_marked_elements (firedrake/adapt.py) gained PETSc.Log.Event wrappers around each of their previously-unattributed sub-steps (label marking, set_adaptive_parent_label, the adaptLabel transform, Mesh() construction, adaptive_parent_child_cell_maps, re-marking, Netgen re-curving). Without them, -log_view attributed all of this work to refine_marked_elements's own self-time rather than the call that actually spent it, which is what obscured this cost in the first place.

The adaptive_multigrid demo's recorded multigrid iteration counts and convergence plot predated #5215's DMPlexTransform-based mesh adaptivity; both are regenerated here against current refine_marked_elements output.

Tests

tests/firedrake/multigrid/test_adaptive_multigrid.py already exercises adaptive_parent_child_cell_maps across firedrake/netgen meshes and process counts; it passes unchanged at nprocs 1/2/4.

🤖 Generated with Claude Code

… construction

adaptive_parent_child_cell_maps resolved each fine cell's coarse parent
with DMLabelGetValue, called once per fine cell. DMLabelGetValue should
resolve a point in O(1); instead it scans every stratum of the label
linearly. set_adaptive_parent_label gives every coarse cell its own
stratum, so the loop was effectively O(nfine * ncoarse).

Walk the label by stratum instead (DMLabelGetStratumSize/GetStratumIS,
one call per coarse cell), which uses PETSc's O(1) value->stratum hash
map and touches each fine cell exactly once: O(nfine + ncoarse).

Also wrap refine_marked_elements's previously-unattributed sub-steps in
PETSc.Log.Event so -log_view attributes time to the call that spent it
instead of lumping it into refine_marked_elements's own self-time.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@pbrubeck pbrubeck added the base:main Run this PR using a main (dev) build label Aug 5, 2026
Comment thread firedrake/adapt.py Outdated
pbrubeck and others added 3 commits August 5, 2026 15:31
The recorded multigrid iteration counts and convergence plot predate
#5215's DMPlexTransform-based mesh adaptivity; regenerate both against
current refine_marked_elements output.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Describe why walking by stratum is fast rather than reciting the
removed point-wise loop's complexity.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@pbrubeck
pbrubeck requested a review from connorjward August 5, 2026 15:15
Comment thread demos/adaptive_multigrid/adaptive_convergence.png
@pbrubeck

pbrubeck commented Aug 6, 2026

Copy link
Copy Markdown
Contributor Author

Thanks. CI failures are known. Merging

@pbrubeck
pbrubeck merged commit 4ae52cc into main Aug 6, 2026
6 of 7 checks passed
@pbrubeck
pbrubeck deleted the pbrubeck/adaptive-optimise branch August 6, 2026 09:01
leo-collins pushed a commit that referenced this pull request Aug 6, 2026
… construction (#5326)

* Fix quadratic cost in adaptive mesh refinement's coarse/fine cell-map construction
haldaas added a commit to haldaas/firedrake that referenced this pull request Aug 13, 2026
Submesh(mesh, comm=COMM_SELF) cost O(P^2) in the number of mesh points P.
On three ranks a UnitSquareMesh(181, 181) needed 8.9 seconds to build one
serial submesh. Four times more cells cost about ten times more time.

Every preconditioner that builds a subdomain matrix through
firedrake.preconditioners.matis pays this cost, because local_mesh builds
one serial submesh per process.

The cost was in submesh_correct_entity_classes. A submesh on COMM_SELF has
no neighbour, so it has no ghost points and no owned points, and every
point is core. The function set that state one point at a time.
DMLabelSetValue moves the stratum into a hash set and destroys its index
set. The DMLabelHasPoint of the next iteration reads the hash set back,
sorts it, and rebuilds the index set and the bit array. Each point
therefore cost O(k log k) in the size k of the core stratum so far.

Set the three strata in bulk instead. This is the same fix as the one in
PR firedrakeproject#5326 for the coarse and fine cell maps.

submesh_create also built a temporary label that marked every cell, and
passed it to DMPlexFilter. DMPlexFilter selects every cell by itself when
it gets no label, so the label was unnecessary. Submesh now leaves
label_name as None for a codim-0 submesh of the whole mesh. A codim-1
submesh still passes the "depth" label, because DMPlexFilter with no label
selects cells, not facets.

Measured on three ranks, for submesh_correct_entity_classes alone:

    cells   before     after
      717    34 ms   0.12 ms
     2787   320 ms   0.43 ms
    11059  3214 ms   1.85 ms
    22023  8907 ms   3.80 ms

The whole submesh build at 11059 cells falls from 3256 ms to 71 ms.

Both paths produce the same submesh. On UnitSquareMesh(128, 128) with
three ranks, the submesh built with a label and the submesh built without
one have the same chart, the same subpoint index set, and the same
pyop2_core, pyop2_owned, pyop2_ghost and Face Sets strata.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pbrubeck added a commit that referenced this pull request Aug 18, 2026
… construction (#5326)

* Fix quadratic cost in adaptive mesh refinement's coarse/fine cell-map construction
pbrubeck added a commit that referenced this pull request Aug 19, 2026
* Make serial submesh construction linear in the number of points

Submesh(mesh, comm=COMM_SELF) cost O(P^2) in the number of mesh points P.
On three ranks a UnitSquareMesh(181, 181) needed 8.9 seconds to build one
serial submesh. Four times more cells cost about ten times more time.

Every preconditioner that builds a subdomain matrix through
firedrake.preconditioners.matis pays this cost, because local_mesh builds
one serial submesh per process.

The cost was in submesh_correct_entity_classes. A submesh on COMM_SELF has
no neighbour, so it has no ghost points and no owned points, and every
point is core. The function set that state one point at a time.
DMLabelSetValue moves the stratum into a hash set and destroys its index
set. The DMLabelHasPoint of the next iteration reads the hash set back,
sorts it, and rebuilds the index set and the bit array. Each point
therefore cost O(k log k) in the size k of the core stratum so far.

Set the three strata in bulk instead. This is the same fix as the one in
PR #5326 for the coarse and fine cell maps.

submesh_create also built a temporary label that marked every cell, and
passed it to DMPlexFilter. DMPlexFilter selects every cell by itself when
it gets no label, so the label was unnecessary. Submesh now leaves
label_name as None for a codim-0 submesh of the whole mesh. A codim-1
submesh still passes the "depth" label, because DMPlexFilter with no label
selects cells, not facets.

Measured on three ranks, for submesh_correct_entity_classes alone:

    cells   before     after
      717    34 ms   0.12 ms
     2787   320 ms   0.43 ms
    11059  3214 ms   1.85 ms
    22023  8907 ms   3.80 ms

The whole submesh build at 11059 cells falls from 3256 ms to 71 ms.

Both paths produce the same submesh. On UnitSquareMesh(128, 128) with
three ranks, the submesh built with a label and the submesh built without
one have the same chart, the same subpoint index set, and the same
pyop2_core, pyop2_owned, pyop2_ghost and Face Sets strata.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

base:main Run this PR using a main (dev) build

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants