From d647ae319b153804263d3fca9ebbc583c12d9bf3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 4 Aug 2026 13:17:25 +0100 Subject: [PATCH] Fix coarse-to-fine cell maps on an adaptively refined coarse mesh coarse_to_fine_cells assumed that the cells kept by the non-overlapped plex were the leading cells of the overlapped one. That holds for a plex distributed by DMPlexDistributeOverlap, which appends its halo, but not for one built by DMPlexTransform, which numbers cells by refinement case and so interleaves owned with halo cells. Mask the absent cells and scatter each surviving entry to its non-overlapped number, rather than dropping them and slicing. Rewrite the docstring and comments of coarse_to_fine_cells to name the three cell numberings it moves between. Co-Authored-By: Claude Opus 5 --- firedrake/cython/mgimpl.pyx | 86 +++++++++++++------ .../multigrid/test_adaptive_multigrid.py | 32 +++++++ 2 files changed, 93 insertions(+), 25 deletions(-) diff --git a/firedrake/cython/mgimpl.pyx b/firedrake/cython/mgimpl.pyx index ac250959e8..4631c32fcd 100644 --- a/firedrake/cython/mgimpl.pyx +++ b/firedrake/cython/mgimpl.pyx @@ -381,14 +381,44 @@ def adaptive_parent_child_cell_maps(PETSc.DM coarse_dm, @cython.boundscheck(False) @cython.wraparound(False) def coarse_to_fine_cells(mc, mf, clgmaps, flgmaps): - """Return a map from (renumbered) cells in a coarse mesh to those - in a refined fine mesh. - - :arg mc: the coarse mesh to create the map from. - :arg mf: the fine mesh to map to. - :arg clgmaps: coarse lgmaps (non-overlapped and overlapped) - :arg flgmaps: fine lgmaps (non-overlapped and overlapped) - :returns: Two arrays, one mapping coarse to fine cells, the second fine to coarse cells. + """Map the cells of a coarse mesh to those of its uniform refinement. + + Parameters + ---------- + mc : MeshGeometry + The coarse mesh. + mf : MeshGeometry + The fine mesh, obtained by uniformly refining the non-overlapped + plex of ``mc``. + clgmaps : tuple + The coarse ``(non-overlapped, overlapped)`` point local-to-global maps. + flgmaps : tuple + The fine ``(non-overlapped, overlapped)`` point local-to-global maps. + + Returns + ------- + numpy.ndarray + Map from each owned coarse cell to the fine cells it was split into. + numpy.ndarray + Map from each owned fine cell to the coarse cell it came from. + + Notes + ----- + Three numberings of the same cells meet here: + + 1. Firedrake numbering, which lists owned cells before halo cells. The + returned maps are indexed by, and contain, these numbers. + 2. Overlapped plex numbering, that of ``mesh.topology_dm``. + `get_entity_renumbering` translates between 1 and 2. + 3. Non-overlapped plex numbering, that of the halo-free plex that was + refined. Only here does the parent relation hold: uniform refinement + splits cell ``p`` into cells ``p*nref`` to ``p*nref + nref - 1``. + + Applying an overlapped local-to-global map and then a non-overlapped + global-to-local one translates between 2 and 3. Those two numberings are + genuinely different orders, not a common prefix plus a halo: a plex built + by ``DMPlexTransform`` (an adaptively refined one) numbers its cells by + refinement case, so its owned cells are interleaved with its halo cells. """ cdef: PETSc.DM cdm, fdm @@ -396,7 +426,7 @@ def coarse_to_fine_cells(mc, mf, clgmaps, flgmaps): PetscInt i, ccell, fcell, nfine np.ndarray coarse_to_fine np.ndarray fine_to_coarse - np.ndarray co2n, fn2o, idx + np.ndarray co2n, fn2o, idx, found, permuted cdm = mc.topology_dm fdm = mf.topology_dm @@ -404,6 +434,8 @@ def coarse_to_fine_cells(mc, mf, clgmaps, flgmaps): nref = 2 ** dim ncoarse = mc.cell_set.size nfine = mf.cell_set.size + # co2n: coarse overlapped plex cell -> coarse Firedrake cell + # fn2o: fine Firedrake cell -> fine overlapped plex cell co2n, _ = get_entity_renumbering(cdm, mc._cell_numbering, "cell") _, fn2o = get_entity_renumbering(fdm, mf._cell_numbering, "cell") coarse_to_fine = np.full((ncoarse, nref), -1, dtype=PETSc.IntType) @@ -411,34 +443,38 @@ def coarse_to_fine_cells(mc, mf, clgmaps, flgmaps): # Walk owned fine cells: cStart, cEnd = 0, nfine + # In serial the overlapped and non-overlapped plexes are the same plex, + # so both maps already speak the numbering the parent relation holds in. if mc.comm.size > 1: + # Cells are the leading points of a plex chart, so these point maps + # can be applied to cell numbers directly. cno, co = clgmaps fno, fo = flgmaps - # Compute global numbers of original cell numbers + # Rebase fn2o onto the fine non-overlapped plex, one map per arrow: + # fine Firedrake cell -> overlapped -> global -> non-overlapped. fo.apply(fn2o, result=fn2o) - # Compute local numbers of original cells on non-overlapped mesh fn2o = fno.applyInverse(fn2o, PETSc.LGMap.MapMode.MASK) - # Need to permute order of co2n so it maps from non-overlapped - # cells to new cells (these may have changed order). Need to - # map all known cells through. + # Rebase co2n the same way, but here it is the *index* that changes + # numbering, not the value, so send every local coarse cell through + # the translation. MASK gives -1 for cells the non-overlapped plex + # does not have. idx = np.arange(mc.cell_set.total_size, dtype=PETSc.IntType) - # LocalToGlobal co.apply(idx, result=idx) - # GlobalToLocal - # Drop values that did not exist on non-overlapped mesh - idx = cno.applyInverse(idx, PETSc.LGMap.MapMode.DROP) - co2n = co2n[idx] + idx = cno.applyInverse(idx, PETSc.LGMap.MapMode.MASK) + # idx[i] is where overlapped cell i lands, so scatter rather than + # slice: the surviving cells need not be the leading ones. + found = idx >= 0 + permuted = np.empty(found.sum(), dtype=PETSc.IntType) + permuted[idx[found]] = co2n[found] + co2n = permuted for c in range(cStart, cEnd): - # get original (overlapped) cell number + # Every owned fine cell exists on the non-overlapped plex. fcell = fn2o[c] - # The owned cells should map into non-overlapped cell numbers - # (due to parallel growth strategy) assert 0 <= fcell < cEnd - # Find original coarse cell (fcell / nref) and then map - # forward to renumbered coarse cell (again non-overlapped - # cells should map into owned coarse cells) + # Uniform refinement numbers the nref children of a cell + # consecutively, so integer division recovers the parent. ccell = co2n[fcell // nref] assert 0 <= ccell < ncoarse fine_to_coarse[c, 0] = ccell diff --git a/tests/firedrake/multigrid/test_adaptive_multigrid.py b/tests/firedrake/multigrid/test_adaptive_multigrid.py index e333c3be74..557b18889f 100644 --- a/tests/firedrake/multigrid/test_adaptive_multigrid.py +++ b/tests/firedrake/multigrid/test_adaptive_multigrid.py @@ -297,6 +297,38 @@ def test_adapt_after_uniform_refinement(coarse_mesh, refine): _assert_adapt_after_uniform_refinement(mh) +@pytest.mark.parallel([1, 2, 4]) +@pytest.mark.parametrize("refine", [1, 2]) +def test_adapt_before_uniform_refinement(coarse_mesh, refine): + """An adaptively refined mesh can be uniformly refined into a hierarchy. + Its plex numbers cells by refinement case, so its owned cells are + interleaved with its halo cells, which the cell maps must not assume away. + """ + netgen_flags = {} if hasattr(coarse_mesh, "netgen_mesh") else None + + M = FunctionSpace(coarse_mesh, "DG", 0) + markers = Function(M) + markers.dat.data_wo[:1] = 1 + mesh = coarse_mesh.refine_marked_elements(markers) + + mh = MeshHierarchy(mesh, refine, netgen_flags=netgen_flags) + assert len(mh) == refine + 1 + assert np.allclose(assemble(1*dx(mh[-1])), assemble(1*dx(coarse_mesh))) + + nref = 2 ** mesh.topological_dimension + for level in range(refine): + coarse_to_fine = mh.coarse_to_fine_cells[level] + fine_to_coarse = mh.fine_to_coarse_cells[level + 1] + assert coarse_to_fine.shape == (mh[level].cell_set.size, nref) + assert fine_to_coarse.shape == (mh[level + 1].cell_set.size, 1) + # Uniform refinement splits every owned coarse cell into nref owned + # fine cells, each of which points back at the cell it came from. + assert (coarse_to_fine >= 0).all() + assert (fine_to_coarse >= 0).all() + parents = np.arange(coarse_to_fine.shape[0]).reshape(-1, 1) + assert (fine_to_coarse[coarse_to_fine, 0] == parents).all() + + @pytest.mark.parallel([1, 2, 4]) @pytest.mark.parametrize("operator", ["prolong", "inject"]) def test_DG0(mh, operator):