diff --git a/demos/adaptive_multigrid/adaptive_convergence.png b/demos/adaptive_multigrid/adaptive_convergence.png index 295d8de024..b1fe43cde2 100644 Binary files a/demos/adaptive_multigrid/adaptive_convergence.png and b/demos/adaptive_multigrid/adaptive_convergence.png differ diff --git a/demos/adaptive_multigrid/adaptive_multigrid.py.rst b/demos/adaptive_multigrid/adaptive_multigrid.py.rst index a8ed630a08..81276fcf44 100644 --- a/demos/adaptive_multigrid/adaptive_multigrid.py.rst +++ b/demos/adaptive_multigrid/adaptive_multigrid.py.rst @@ -235,18 +235,18 @@ Moreover, the multigrid iteration count is robust to the level of refinement :: 0 2 1 8 2 8 - 3 8 - 4 8 - 5 8 - 6 8 - 7 8 - 8 8 - 9 9 - 10 9 - 11 9 - 12 9 - 13 9 - 14 9 + 3 7 + 4 7 + 5 7 + 6 7 + 7 7 + 8 7 + 9 7 + 10 7 + 11 7 + 12 7 + 13 7 + 14 7 ======== ================ A runnable python version of this demo can be found :demo:`here`. diff --git a/firedrake/adapt.py b/firedrake/adapt.py index 071ca8231a..f1b6a4851b 100644 --- a/firedrake/adapt.py +++ b/firedrake/adapt.py @@ -9,6 +9,7 @@ from firedrake.functionspace import FunctionSpace from firedrake.mesh import Mesh, DISTRIBUTION_PARAMETERS_NOOP from firedrake.netgen import _transfer_high_order_coordinates +from firedrake.petsc import PETSc # PETSc's DMAdaptFlag value requesting refinement, for the adapt label. @@ -28,20 +29,22 @@ def _adapt_marked_cells(mesh, cell_marker): dm = mesh.topology_dm ncoarse = mesh.cell_set.size - dm.createLabel(ADAPT_LABEL) - adapt_label = dm.getLabel(ADAPT_LABEL) - adapt_indicator = np.zeros(cell_marker.dat.data_ro_with_halos.shape, dtype=IntType) - adapt_indicator[:ncoarse] = cell_marker.dat.data_ro.real > 0 - dmcommon.mark_points_with_function_array( - dm, cell_marker.function_space().dm.getSection(), 0, - adapt_indicator, adapt_label, DM_ADAPT_REFINE, - ) + with PETSc.Log.Event("AdaptiveRefine: mark cells"): + dm.createLabel(ADAPT_LABEL) + adapt_label = dm.getLabel(ADAPT_LABEL) + adapt_indicator = np.zeros(cell_marker.dat.data_ro_with_halos.shape, dtype=IntType) + adapt_indicator[:ncoarse] = cell_marker.dat.data_ro.real > 0 + dmcommon.mark_points_with_function_array( + dm, cell_marker.function_space().dm.getSection(), 0, + adapt_indicator, adapt_label, DM_ADAPT_REFINE, + ) parameters = {"dm_plex_transform_type": "refine_sbr"} try: # options_prefix="" is essential with petsctools.inserted_options(parameters=parameters, options_prefix=""): - new_dm = dm.adaptLabel(ADAPT_LABEL) + with PETSc.Log.Event("AdaptiveRefine: adaptLabel"): + new_dm = dm.adaptLabel(ADAPT_LABEL) finally: # Ensure the temporary label is removed even if adaptation fails dm.removeLabel(ADAPT_LABEL) @@ -99,33 +102,37 @@ def refine_marked_elements(mesh, cell_marker): num_refinements = max(int(np.rint(num_refinements)), 1) coarse_dm = mesh.topology_dm - impl.set_adaptive_parent_label(coarse_dm, mesh._cell_numbering, PARENT_LABEL) + with PETSc.Log.Event("AdaptiveRefine: set_adaptive_parent_label"): + impl.set_adaptive_parent_label(coarse_dm, mesh._cell_numbering, PARENT_LABEL) current_mesh = mesh current_mark = cell_marker try: for ref in range(num_refinements): new_dm = _adapt_marked_cells(current_mesh, current_mark) - current_mesh = Mesh( - new_dm, - dim=mesh.geometric_dimension, - reorder=False, - distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, - comm=mesh.comm, - tolerance=mesh.tolerance, - ) - coarse_to_fine, fine_to_coarse = impl.adaptive_parent_child_cell_maps( - coarse_dm, new_dm, current_mesh._cell_numbering, PARENT_LABEL - ) + with PETSc.Log.Event("AdaptiveRefine: Mesh()"): + current_mesh = Mesh( + new_dm, + dim=mesh.geometric_dimension, + reorder=False, + distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, + comm=mesh.comm, + tolerance=mesh.tolerance, + ) + with PETSc.Log.Event("AdaptiveRefine: adaptive_parent_child_cell_maps"): + coarse_to_fine, fine_to_coarse = impl.adaptive_parent_child_cell_maps( + coarse_dm, new_dm, current_mesh._cell_numbering, PARENT_LABEL + ) if ref < num_refinements - 1: - # A cell asking for n refinements stays marked until n rounds - # have happened, so its descendants inherit n minus the number - # of rounds so far. - ancestor = fine_to_coarse[:, 0] - refined = ancestor >= 0 - current_mark = Function(FunctionSpace(current_mesh, "DG", 0)) - current_mark.dat.data_wo[refined] = \ - cell_marker.dat.data_ro[ancestor[refined]] - (ref + 1) + with PETSc.Log.Event("AdaptiveRefine: re-mark"): + # A cell asking for n refinements stays marked until n rounds + # have happened, so its descendants inherit n minus the number + # of rounds so far. + ancestor = fine_to_coarse[:, 0] + refined = ancestor >= 0 + current_mark = Function(FunctionSpace(current_mesh, "DG", 0)) + current_mark.dat.data_wo[refined] = \ + cell_marker.dat.data_ro[ancestor[refined]] - (ref + 1) finally: # Ensure the temporary label is removed even if adaptation fails coarse_dm.removeLabel(PARENT_LABEL) @@ -134,7 +141,8 @@ def refine_marked_elements(mesh, cell_marker): if hasattr(mesh, "netgen_mesh"): order = mesh.coordinates.function_space().ufl_element().degree() if order > 1: - final_mesh = _transfer_high_order_coordinates(mesh, final_mesh, order) + with PETSc.Log.Event("AdaptiveRefine: recurve netgen coords"): + final_mesh = _transfer_high_order_coordinates(mesh, final_mesh, order) final_mesh.topology_dm.removeLabel(PARENT_LABEL) final_mesh.adaptive_parent = mesh diff --git a/firedrake/cython/mgimpl.pyx b/firedrake/cython/mgimpl.pyx index 4631c32fcd..4b3a9a3e77 100644 --- a/firedrake/cython/mgimpl.pyx +++ b/firedrake/cython/mgimpl.pyx @@ -296,8 +296,10 @@ def adaptive_parent_child_cell_maps(PETSc.DM coarse_dm, cdef: PetscInt ncoarse = num_owned_cells(coarse_dm) PetscInt nfine = num_owned_cells(fine_dm) - PetscInt cStart, cEnd, c, off, parent, max_children + PetscInt cStart, cEnd, c, off, parent, i, stratum_size, max_children DMLabel parent_label = NULL + PETSc.PetscIS stratum_is = NULL + const PetscInt *stratum_points = NULL PetscInt[::1] child_counts PetscInt[:, ::1] coarse_to_fine PetscInt[:, ::1] fine_to_coarse @@ -307,14 +309,26 @@ def adaptive_parent_child_cell_maps(PETSc.DM coarse_dm, fine_to_coarse = np.full((nfine, 1), -1, dtype=IntType) child_counts = np.zeros(ncoarse, dtype=IntType) cStart, cEnd = fine_dm.getHeightStratum(0) - for c in range(cStart, cEnd): - CHKERR(PetscSectionGetOffset(fine_cell_numbering.sec, c, &off)) - if not (0 <= off < nfine): + # Walking by stratum (coarse cell) resolves each one through PETSc's O(1) + # value -> stratum hash map and touches every fine cell exactly once, for + # O(nfine + ncoarse) overall. + for parent in range(ncoarse): + CHKERR(DMLabelGetStratumSize(parent_label, parent, &stratum_size)) + if stratum_size <= 0: continue - CHKERR(DMLabelGetValue(parent_label, c, &parent)) - if 0 <= parent < ncoarse: + CHKERR(DMLabelGetStratumIS(parent_label, parent, &stratum_is)) + CHKERR(ISGetIndices(stratum_is, &stratum_points)) + for i in range(stratum_size): + c = stratum_points[i] + if not (cStart <= c < cEnd): + continue + CHKERR(PetscSectionGetOffset(fine_cell_numbering.sec, c, &off)) + if not (0 <= off < nfine): + continue fine_to_coarse[off, 0] = parent child_counts[parent] += 1 + CHKERR(ISRestoreIndices(stratum_is, &stratum_points)) + CHKERR(ISDestroy(&stratum_is)) # coarse_to_fine is rectangular, so every coarse cell's row must be wide # enough for its most prolific sibling. Different coarse cells can be