diff --git a/AGENTS.md b/AGENTS.md index de62cc1d29..ffc2ca2295 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -74,8 +74,21 @@ toolchain: otherwise touching the same whole array repeatedly outside of a single vectorized expression — that is what defeats NumPy's own performance model, on top of bypassing PyOP2/code-generation for mesh-bound data. -* **Docstrings:** All public-facing APIs must include properly formatted `numpydoc`-style docstrings. +* **Docstrings Are Always `numpydoc`:** Every docstring you write or touch — public API, private helper, + Cython function in `firedrake/cython/*.pyx`, test helper — must be `numpydoc`, using its section + headings (`Parameters`, `Returns`, `Raises`, `Notes`). Never write the old Sphinx field-list style + (`:arg x:`, `:param x:`, `:returns:`, `:rtype:`) in new or edited code, and do not copy it from the + surrounding file: much of Firedrake predates the convention, so matching the neighbouring docstrings + is precisely the wrong instinct — this is the one place where "preserve the existing style" does not + apply. Being private, internal, or compiled is not an excuse to skip the docstring, to downgrade its + format, or to leave the arguments undocumented: give every parameter and every return value its + `numpydoc` entry, however small the helper. * **Type Hints:** New code should include type hints on function/method signatures. +* **Demos Are Literate Programs:** `pylit` converts each `demos//.py.rst` into a `.py` that + `tests/firedrake/demos/test_demos_run.py` executes, so prose and code must stay in step. A paragraph + ending in `::` makes the indented block after it *executable*; a `.. code-block:: python` directive is + excluded from that rule, so its snippet renders in the docs but never runs. Prefer `::` — reach for + the directive only for an illustrative fragment naming things the demo never defines. ## Testing Requirements diff --git a/demos/adaptive_multigrid/adaptive_multigrid.py.rst b/demos/adaptive_multigrid/adaptive_multigrid.py.rst index 44ff4eac0f..a8ed630a08 100644 --- a/demos/adaptive_multigrid/adaptive_multigrid.py.rst +++ b/demos/adaptive_multigrid/adaptive_multigrid.py.rst @@ -1,13 +1,12 @@ -Adaptive Multigrid Methods using AdaptiveMeshHierarchy -====================================================== +Adaptive Multigrid Methods +========================== Contributed by Anurag Rao. The purpose of this demo is to show how to use Firedrake's multigrid solver on a hierarchy of adaptively refined Netgen meshes. -We will first have a look at how to use the :class:`~.AdaptiveMeshHierarchy` to construct the mesh hierarchy with Netgen meshes, then we will consider a solution to the Poisson problem on an L-shaped domain. -Finally, we will show how to use the :class:`~.AdaptiveMeshHierarchy` and :class:`~.AdaptiveTransferManager` to construct a scalable solver. The :class:`~.AdaptiveMeshHierarchy` contains information of the mesh hierarchy and the parent child relations between the meshes. -The :class:`~.AdaptiveTransferManager` deals with the transfer operator logic across any given levels in the hierarchy. +A :func:`~.MeshHierarchy` is not restricted to uniform refinement: the same object records the parent child relations between adaptively refined meshes, and grows a level at a time as the solution is resolved. +We will first have a look at how to construct such a hierarchy from Netgen meshes, then we will consider a solution to the Poisson problem on an L-shaped domain, and finally we will use the hierarchy to construct a scalable solver. We begin by importing the necessary libraries :: from firedrake import * @@ -28,16 +27,15 @@ We begin with the L-shaped domain, which we build as the union of two rectangles ngmsh = geo.GenerateMesh(maxh=0.5) mesh = Mesh(ngmsh) -It is important to convert the initial Netgen mesh into a Firedrake mesh before constructing the :class:`~.AdaptiveMeshHierarchy`. To call the constructor to the hierarchy, we must pass the initial mesh. Our initial mesh looks like this: +It is important to convert the initial Netgen mesh into a Firedrake mesh before constructing the :func:`~.MeshHierarchy`. To call the constructor to the hierarchy, we must pass the initial mesh. Our initial mesh looks like this: .. figure:: initial_mesh.png :align: center :alt: Initial mesh. -We will also initialize the :class:`~.AdaptiveTransferManager` here: :: +We initialize the :func:`~.MeshHierarchy` here. The default of zero uniform refinement levels gives a hierarchy holding just the initial mesh, which we will grow adaptively below; passing a positive number instead would start us off with that many uniformly refined levels, and the adaptive levels would stack on top of them just the same: :: - amh = AdaptiveMeshHierarchy(mesh) - atm = AdaptiveTransferManager() + mh = MeshHierarchy(mesh) Poisson Problem --------------- @@ -62,15 +60,12 @@ Our approach strongly follows the similar problem in this `lecture course 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) + finally: + # Ensure the temporary label is removed even if adaptation fails + dm.removeLabel(ADAPT_LABEL) + + # The transform propagates every label, including the temporary adapt + # label and the coarse mesh's stale pyop2_core/owned/ghost point + # classification. Mesh() skips recomputing that classification if it's + # already present, so it must be dropped here to force a fresh one for + # the new mesh's own point count and distribution. + for label in ("pyop2_core", "pyop2_owned", "pyop2_ghost", ADAPT_LABEL): + if new_dm.hasLabel(label): + new_dm.removeLabel(label) + + return new_dm + + +def _copy_adaptive_refinement_metadata(source_mesh, target_mesh): + """Copy mesh-construction metadata from a mesh onto its adaptively-derived successor.""" + target_mesh._distribution_parameters = dict(source_mesh._distribution_parameters) + target_mesh._did_reordering = source_mesh._did_reordering + target_mesh._tolerance = source_mesh.tolerance + if hasattr(source_mesh, "netgen_mesh") and not hasattr(target_mesh, "netgen_mesh"): + target_mesh.netgen_mesh = source_mesh.netgen_mesh + if hasattr(source_mesh, "netgen_flags") and not hasattr(target_mesh, "netgen_flags"): + target_mesh.netgen_flags = source_mesh.netgen_flags + + +def refine_marked_elements(mesh, cell_marker): + """Adaptively refine a mesh using a DG0 marking function. + + Positive integer marker values request repeated refinement of the + corresponding cells. Curved Netgen meshes are re-curved to the + original coordinate degree after refinement. + + Parameters + ---------- + mesh + The mesh to refine. + cell_marker + A DG0 `~firedrake.function.Function` on ``mesh``: cells with a + positive value ``n`` are refined ``n`` times. + + Returns + ------- + MeshGeometry + The adaptively refined mesh, with ``adaptive_parent`` set to + ``mesh`` and ``adaptive_cell_maps`` set to the + ``(coarse_to_fine, fine_to_coarse)`` cell maps relative to it. + + """ + with cell_marker.dat.vec_ro as v: + _, num_refinements = v.max() + # Always run at least one adaptation pass, even when no cell is marked, + # so that a fresh mesh (with its own cell maps) is produced uniformly. + 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) + + 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 + ) + 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) + finally: + # Ensure the temporary label is removed even if adaptation fails + coarse_dm.removeLabel(PARENT_LABEL) + + final_mesh = current_mesh + 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) + + final_mesh.topology_dm.removeLabel(PARENT_LABEL) + final_mesh.adaptive_parent = mesh + final_mesh.adaptive_cell_maps = (coarse_to_fine, fine_to_coarse) + _copy_adaptive_refinement_metadata(mesh, final_mesh) + return final_mesh diff --git a/firedrake/cython/mgimpl.pyx b/firedrake/cython/mgimpl.pyx index b9b41bd32f..ac250959e8 100644 --- a/firedrake/cython/mgimpl.pyx +++ b/firedrake/cython/mgimpl.pyx @@ -6,6 +6,7 @@ import numpy as np from firedrake.cython import dmcommon from firedrake.petsc import PETSc from firedrake.utils import IntType +from pyop2.mpi import MPI cimport numpy as np cimport petsc4py.PETSc as PETSc @@ -97,6 +98,9 @@ def coarse_to_fine_nodes(Vc, Vf, np.ndarray coarse_to_fine_cells): k = 0 for l in range(fine_cell_per_coarse_cell): fine = coarse_to_fine_cells[i, l] + if fine < 0: + k += fine_per_cell * ratio + continue for layer in range(ratio): fine_layer = coarse_layer * ratio + layer for m in range(fine_per_cell): @@ -107,6 +111,9 @@ def coarse_to_fine_nodes(Vc, Vf, np.ndarray coarse_to_fine_cells): k = 0 for l in range(fine_cell_per_coarse_cell): fine = coarse_to_fine_cells[i, l] + if fine < 0: + k += fine_per_cell + continue for m in range(fine_per_cell): coarse_to_fine_map[node, k] = fine_map[fine, m] k += 1 @@ -149,6 +156,8 @@ def fine_to_coarse_nodes(Vf, Vc, np.ndarray fine_to_coarse_cells): for i in range(fine_cells): for l, coarse_cell in enumerate(fine_to_coarse_cells[i, :]): + if coarse_cell < 0: + continue for j in range(fine_per_cell): node = fine_map[i, j] if extruded: @@ -192,6 +201,144 @@ def create_lgmap(PETSc.DM dm): return lgmap +cdef PetscInt num_owned_cells(PETSc.DM dm) except? -1: + """Number of cells this rank owns, i.e. the number of Firedrake cell + numbers the DM's cell numbering hands out to non-ghost cells. + + Parameters + ---------- + dm : PETSc.DM + The DMPlex encapsulating the mesh topology, with its PyOP2 entity + classes already marked. + + Returns + ------- + PetscInt + The number of core plus owned cells. + + """ + return dmcommon.get_entity_classes(dm)[dm.getDimension(), 1] + + +@cython.boundscheck(False) +@cython.wraparound(False) +def set_adaptive_parent_label(PETSc.DM coarse_dm, + PETSc.Section coarse_cell_numbering, + label_name): + """Seed each coarse cell's own Firedrake cell number onto a DMPlex label. + + Must be called *before* refining ``coarse_dm``. Since the refinement + transform propagates labels from a cell to its children, every cell of + every subsequent refinement then carries the number of the coarse cell it + descends from, which `adaptive_parent_child_cell_maps` reads back. + + Parameters + ---------- + coarse_dm : PETSc.DM + The coarse, pre-refinement, mesh DMPlex. + coarse_cell_numbering : PETSc.Section + The coarse mesh's cell numbering section. + label_name : str + Name of the label to create on ``coarse_dm`` and populate with each + owned cell's Firedrake cell number. Any existing label of that name + is discarded. + + """ + cdef: + PetscInt ncoarse = num_owned_cells(coarse_dm) + PetscInt cStart, cEnd, c, off + DMLabel parent_label = NULL + + if coarse_dm.hasLabel(label_name): + coarse_dm.removeLabel(label_name) + coarse_dm.createLabel(label_name) + label_name = label_name.encode() + CHKERR(DMGetLabel(coarse_dm.dm, label_name, &parent_label)) + cStart, cEnd = coarse_dm.getHeightStratum(0) + for c in range(cStart, cEnd): + CHKERR(PetscSectionGetOffset(coarse_cell_numbering.sec, c, &off)) + if 0 <= off < ncoarse: + CHKERR(DMLabelSetValue(parent_label, c, off)) + + +@cython.boundscheck(False) +@cython.wraparound(False) +def adaptive_parent_child_cell_maps(PETSc.DM coarse_dm, + PETSc.DM fine_dm, + PETSc.Section fine_cell_numbering, + label_name): + """Build Firedrake-numbered parent/child cell maps from a DMPlex label. + + ``fine_dm`` must be a DMPlex obtained by refining, however many times, the + ``coarse_dm`` that `set_adaptive_parent_label` was seeded on. + + Parameters + ---------- + coarse_dm : PETSc.DM + The coarse, parent, mesh DMPlex. + fine_dm : PETSc.DM + The refined, child, mesh DMPlex. + fine_cell_numbering : PETSc.Section + The fine mesh's cell numbering section. + label_name : str + Name of the label on ``fine_dm``, propagated from + `set_adaptive_parent_label`, mapping each fine cell to its coarse + parent's Firedrake cell number. + + Returns + ------- + tuple of numpy.ndarray + The ``(coarse_to_fine, fine_to_coarse)`` Firedrake-numbered cell maps, + padded with -1 where a coarse cell has fewer children than the + busiest one. + + """ + cdef: + PetscInt ncoarse = num_owned_cells(coarse_dm) + PetscInt nfine = num_owned_cells(fine_dm) + PetscInt cStart, cEnd, c, off, parent, max_children + DMLabel parent_label = NULL + PetscInt[::1] child_counts + PetscInt[:, ::1] coarse_to_fine + PetscInt[:, ::1] fine_to_coarse + + label_name = label_name.encode() + CHKERR(DMGetLabel(fine_dm.dm, label_name, &parent_label)) + 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): + continue + CHKERR(DMLabelGetValue(parent_label, c, &parent)) + if 0 <= parent < ncoarse: + fine_to_coarse[off, 0] = parent + child_counts[parent] += 1 + + # 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 + # refined a different number of times, so this varies by process; take + # the max across all ranks so the array shape agrees everywhere. + max_children = 0 + for c in range(ncoarse): + if child_counts[c] > max_children: + max_children = child_counts[c] + max_children = fine_dm.comm.tompi4py().allreduce(max_children, op=MPI.MAX) + coarse_to_fine = np.full((ncoarse, max_children), -1, dtype=IntType) + # Re-walk the fine cells (in Firedrake order this time, via + # fine_to_coarse) appending each one to its parent's row. child_counts is + # reused as a per-parent write cursor, reset to zero first. + child_counts[:] = 0 + for c in range(nfine): + parent = fine_to_coarse[c, 0] + if parent >= 0: + coarse_to_fine[parent, child_counts[parent]] = c + child_counts[parent] += 1 + + return np.asarray(coarse_to_fine), np.asarray(fine_to_coarse) + + # Exposition: # # These next functions compute maps from coarse mesh cells to fine diff --git a/firedrake/dmhooks.py b/firedrake/dmhooks.py index 5b14465010..189203e25d 100644 --- a/firedrake/dmhooks.py +++ b/firedrake/dmhooks.py @@ -453,12 +453,12 @@ def _refine_adaptive(dm): Return the DM of the `_SNESContext` reconstructed on the adaptively-refined mesh using `_SNESContext.marking_callback` to mark the cells to be refined. """ - from firedrake.mg.adaptive_hierarchy import AdaptiveMeshHierarchy + from firedrake.mg.mesh import MeshHierarchy from firedrake.mg.ufl_utils import refine from firedrake.mg.utils import get_level # DMAdaptorAdapt() unconditionally destroys its input DM, and each - # adapted input DM remains a level in the AdaptiveMeshHierarchy. + # adapted input DM remains a level in the mesh hierarchy. # Increase the reference count so the coarse DM survives. dm.incRef() @@ -469,11 +469,9 @@ def _refine_adaptive(dm): mesh = current_solution.function_space().mesh() hierarchy, level = get_level(mesh) if hierarchy is None: - hierarchy = AdaptiveMeshHierarchy(mesh) + hierarchy = MeshHierarchy(mesh) level = 0 - if not isinstance(hierarchy, AdaptiveMeshHierarchy): - raise RuntimeError("Adaptive SNES refinement requires an AdaptiveMeshHierarchy") if level+1 != len(hierarchy): raise RuntimeError("Adaptive SNES refinement can only add a mesh on top of the finest level") if ctx._marking_callback is None: diff --git a/firedrake/mesh.py b/firedrake/mesh.py index b2b91d1e25..881cdfe523 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -1432,8 +1432,9 @@ def cell_set(self): size = list(self._entity_classes[self.cell_dimension(), :]) return op2.Set(size, "Cells", comm=self.comm) + @staticmethod @PETSc.Log.EventDecorator() - def _set_partitioner(self, plex, distribute, partitioner_type=None): + def _set_partitioner(plex, distribute, partitioner_type=None): """Set partitioner for (re)distributing underlying plex over comm. :arg distribute: Boolean or (sizes, points)-tuple. If (sizes, point)- @@ -2394,6 +2395,9 @@ def __init__(self, coordinates): self.extruded = isinstance(topology, ExtrudedMeshTopology) self.variable_layers = self.extruded and topology.variable_layers self._base_mesh = None # this is set by extruded meshes in a later step + # these are set by firedrake.adapt.refine_marked_elements + self.adaptive_parent = None + self.adaptive_cell_maps = None self.topology = topology self.geometric_shared_data_cache = defaultdict(dict) @@ -2951,68 +2955,26 @@ def __iter__(self): def unique(self): return self - def refine_marked_elements(self, mark, netgen_flags=None): - """Refine a mesh using a DG0 marking function. - - This method requires that the mesh has been constructed from a - netgen mesh. - - :arg mark: the marking function which is a Firedrake DG0 function - with the number of refinements on each cell. - :arg netgen_flags: the dictionary of flags to be passed to ngsPETSc. + @PETSc.Log.EventDecorator() + def refine_marked_elements(self, mark): + """Adaptively refine a mesh using a DG0 marking function. - It includes the option: - - refine_faces, which is a boolean specifying if you want to refine faces. + Parameters + ---------- + mark + A DG0 `~firedrake.function.Function` on this mesh: cells + with a positive value ``n`` are refined ``n`` times. + Returns + ------- + MeshGeometry + The adaptively refined mesh, recording this mesh as its + ``adaptive_parent`` and the cell maps relative to it as its + ``adaptive_cell_maps``, ready to be passed to + :meth:`~firedrake.mg.mesh.HierarchyBase.add_mesh`. """ - utils.check_netgen_installed() - - if not hasattr(self, "netgen_mesh"): - raise ValueError("Adaptive refinement requires a netgen mesh.") - if netgen_flags is None: - netgen_flags = self.netgen_flags - tdim = self.topological_dimension - if tdim not in {2, 3}: - raise NotImplementedError("No implementation for dimension other than 2 and 3.") - with mark.dat.vec as mvec: - if self.sfBC_orig is None: - cstart, cend = self.topology_dm.getHeightStratum(0) - cellNum = list(map(self._cell_numbering.getOffset, range(cstart, cend))) - mark_np = mvec.getArray()[cellNum] - else: - sfBCInv = self.sfBC_orig.createInverse() - _, mvec0 = self.topology_dm.distributeField(sfBCInv, - self._cell_numbering, - mvec) - mark_np = mvec0.getArray() - max_refs = 0 if mark_np.size == 0 else int(mark_np.max()) - # Create a copy of the netgen mesh - netgen_mesh = self.netgen_mesh.Copy() - refine_faces = netgen_flags.get("refine_faces", False) - for r in range(max_refs): - cells = netgen_mesh.Elements3D() if tdim == 3 else netgen_mesh.Elements2D() - cells.NumPy()["refine"] = (mark_np[:len(cells)] > 0) - if tdim == 3: - faces = netgen_mesh.Elements2D() - faces.NumPy()["refine"] = refine_faces - netgen_mesh.Refine(adaptive=True) - mark_np -= 1 - if r < max_refs - 1: - parents = netgen_mesh.parentelements if tdim == 3 else netgen_mesh.parentsurfaceelements - parents = parents.NumPy()["i"] - num_fine_cells = parents.shape[0] - num_coarse_cells = mark_np.size - indices = np.arange(num_fine_cells, dtype=PETSc.IntType) - while (indices >= num_coarse_cells).any(): - fine_cells = (indices >= num_coarse_cells) - indices[fine_cells] = parents[indices[fine_cells]] - mark_np = mark_np[indices] - - return Mesh(netgen_mesh, - reorder=self._did_reordering, - distribution_parameters=self._distribution_parameters, - comm=self.comm, - netgen_flags=netgen_flags) + from firedrake.adapt import refine_marked_elements + return refine_marked_elements(self, mark) @PETSc.Log.EventDecorator() def curve_field(self, order, permutation_tol=None, cg_field=None): @@ -3485,6 +3447,8 @@ def Mesh(meshfile, **kwargs): comm=mesh.comm) temp.netgen_mesh = mesh.netgen_mesh temp.netgen_flags = mesh.netgen_flags + temp.sfBC = mesh.sfBC + temp.sfBC_orig = mesh.sfBC_orig temp._distribution_parameters = mesh._distribution_parameters temp._did_reordering = mesh._did_reordering mesh = temp diff --git a/firedrake/mg/__init__.py b/firedrake/mg/__init__.py index c73e5c7849..77b65ce6c8 100644 --- a/firedrake/mg/__init__.py +++ b/firedrake/mg/__init__.py @@ -1,12 +1,10 @@ from firedrake.mg.mesh import ( # noqa F401 HierarchyBase, MeshHierarchy, ExtrudedMeshHierarchy, NonNestedHierarchy, SemiCoarsenedExtrudedHierarchy, - SubmeshHierarchy, + SubmeshHierarchy, AdaptiveMeshHierarchy, ) from firedrake.mg.interface import ( # noqa F401 prolong, restrict, inject ) -from firedrake.mg.embedded import TransferManager # noqa F401 +from firedrake.mg.embedded import TransferManager, AdaptiveTransferManager # noqa F401 from firedrake.mg.opencascade_mh import OpenCascadeMeshHierarchy # noqa F401 -from firedrake.mg.adaptive_hierarchy import AdaptiveMeshHierarchy # noqa F401 -from firedrake.mg.adaptive_transfer_manager import AdaptiveTransferManager # noqa: F401 diff --git a/firedrake/mg/adaptive_hierarchy.py b/firedrake/mg/adaptive_hierarchy.py deleted file mode 100644 index fc6c2f728a..0000000000 --- a/firedrake/mg/adaptive_hierarchy.py +++ /dev/null @@ -1,85 +0,0 @@ -from firedrake.mesh import MeshGeometry -from firedrake.cofunction import Cofunction -from firedrake.function import Function -from firedrake.mg import HierarchyBase -from firedrake.mg.utils import set_level - -__all__ = ["AdaptiveMeshHierarchy"] - - -class AdaptiveMeshHierarchy(HierarchyBase): - """ - HierarchyBase for hierarchies of adaptively refined meshes. - - Parameters - ---------- - base_mesh - The coarsest mesh in the hierarchy. - nested: bool - A flag to indicate whether the meshes are nested. - - """ - def __init__(self, base_mesh: MeshGeometry, nested: bool = True): - self.meshes = [] - self._meshes = [] - self.nested = nested - self.add_mesh(base_mesh) - - def add_mesh(self, mesh: MeshGeometry): - """ - Adds a mesh into the hierarchy. - - Parameters - ---------- - mesh - The mesh to be added to the finest level. - """ - level = len(self.meshes) - self._meshes.append(mesh) - self.meshes.append(mesh) - set_level(mesh, self, level) - - def adapt(self, eta: Function | Cofunction, theta: float): - """ - Adds a new mesh to the hierarchy by locally refining the finest mesh - with a simplified variant of Dorfler marking. The finest mesh must - come from a netgen mesh. - - Parameters - ---------- - eta - A DG0 :class:`~firedrake.function.Function` with the local error estimator. - theta - The threshold for marking as a fraction of the maximum error. - - Note - ---- - Dorfler marking involves sorting all of the elements by decreasing - error estimator and taking the minimal set that exceeds some fixed - fraction of the total error. What this code implements is the simpler - variant that doesn't have a proof of convergence (as far as I know) - but works as well in practice. - - """ - if not isinstance(eta, (Function, Cofunction)): - raise TypeError(f"eta must be a Function or Cofunction, not a {type(eta).__name__}") - M = eta.function_space() - if M.finat_element.space_dimension() != 1: - raise ValueError("eta must be a Function or Cofunction in DG0") - mesh = self.meshes[-1] - if M.mesh() is not mesh: - raise ValueError("eta must be defined on the finest mesh of the hierarchy") - - # Take the maximum over all processes - with eta.dat.vec_ro as evec: - _, eta_max = evec.max() - - threshold = theta * eta_max - should_refine = eta.dat.data_ro > threshold - - markers = Function(M) - markers.dat.data_wo[should_refine] = 1 - - refined_mesh = mesh.refine_marked_elements(markers) - self.add_mesh(refined_mesh) - return refined_mesh diff --git a/firedrake/mg/adaptive_transfer_manager.py b/firedrake/mg/adaptive_transfer_manager.py deleted file mode 100644 index 149906b845..0000000000 --- a/firedrake/mg/adaptive_transfer_manager.py +++ /dev/null @@ -1,55 +0,0 @@ -""" -This module contains the AdaptiveTransferManager used to perform -transfer operations on AdaptiveMeshHierarchies -""" -from firedrake.mg.embedded import TransferManager -from firedrake.ufl_expr import action, TrialFunction -from firedrake.interpolation import interpolate - - -__all__ = ("AdaptiveTransferManager",) - - -class AdaptiveTransferManager(TransferManager): - """ - TransferManager for adaptively refined mesh hierarchies - """ - def __init__(self, *, native_transfers=None, use_averaging=True): - if native_transfers is not None: - raise NotImplementedError("Custom transfers not implemented.") - super().__init__(native_transfers=native_transfers, use_averaging=use_averaging) - self.cache = {} - - def get_interpolator(self, Vc, Vf): - from firedrake.assemble import assemble - key = (Vc, Vf) - try: - return self.cache[key] - except KeyError: - Iexpr = interpolate(TrialFunction(Vc), Vf) - # TODO reusable matfree Interpolator - I = assemble(Iexpr, mat_type="aij") - return self.cache.setdefault(key, I) - - def forward(self, uc, uf): - from firedrake.assemble import assemble - Vc = uc.function_space() - Vf = uf.function_space() - I = self.get_interpolator(Vc, Vf) - return assemble(action(I, uc), tensor=uf) - - def adjoint(self, uf, uc): - from firedrake.assemble import assemble - Vc = uc.function_space().dual() - Vf = uf.function_space().dual() - I = self.get_interpolator(Vc, Vf) - return assemble(action(uf, I), tensor=uc) - - def prolong(self, uf, uc): - return self.forward(uf, uc) - - def inject(self, uc, uf): - return self.forward(uc, uf) - - def restrict(self, uc, uf): - return self.adjoint(uc, uf) diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py index 2363a58fe5..9a00c277b1 100644 --- a/firedrake/mg/embedded.py +++ b/firedrake/mg/embedded.py @@ -1,6 +1,7 @@ import firedrake import ufl import finat.ufl +import warnings import weakref from enum import IntEnum from firedrake.petsc import PETSc @@ -8,7 +9,7 @@ from finat.element_factory import create_element from .utils import get_level -__all__ = ("TransferManager", ) +__all__ = ("TransferManager",) class Op(IntEnum): @@ -380,3 +381,12 @@ def restrict(self, source, target): self.DG_inv_mass(VDGt).mult(dgv, dgwork) self.V_DG_mass(Vt, VDGt).multTranspose(dgwork, t) self.cache_dat_versions(Vs_star, Op.RESTRICT, source, target) + + +def AdaptiveTransferManager(*args, **kwargs): + """Deprecated alias for `TransferManager`.""" + warnings.warn( + "The ``AdaptiveTransferManager`` class is deprecated and will be removed in a future release. " + "Please use the ``TransferManager`` class instead.", FutureWarning + ) + return TransferManager(*args, **kwargs) diff --git a/firedrake/mg/interface.py b/firedrake/mg/interface.py index f10234449f..4c544535cf 100644 --- a/firedrake/mg/interface.py +++ b/firedrake/mg/interface.py @@ -230,10 +230,6 @@ def inject(fine, coarse): # Introduce an intermediate quadrature target space Vc = Vc.quadrature_space() - kernel, dg = kernels.inject_kernel(Vf, Vc) - if dg and not hierarchy.nested: - raise NotImplementedError("Sorry, we can't do supermesh projections yet!") - coarsest = coarse.zero() Vcoarsest = coarsest.function_space() meshes = hierarchy._meshes @@ -245,6 +241,9 @@ def inject(fine, coarse): coarse = Function(Vc.reconstruct(mesh=meshes[next_level])) Vc = coarse.function_space() Vf = fine.function_space() + kernel, dg = kernels.inject_kernel(Vf, Vc) + if dg and not hierarchy.nested: + raise NotImplementedError("Multigrid DG injection not implemented on non-nested hierarchies.") if not dg: compose_map = lambda u: utils.coarse_node_to_fine_node_map(Vc, u.function_space()) node_locations = utils.physical_node_locations(Vc) diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py index 749466f05e..70021eec11 100644 --- a/firedrake/mg/kernels.py +++ b/firedrake/mg/kernels.py @@ -352,7 +352,8 @@ def inject_kernel(Vf, Vc): level_ratio = (Vf.mesh().layers - 1) // (Vc.mesh().layers - 1) else: level_ratio = 1 - key = (("inject", level_ratio) + ncandidate = hierarchy.coarse_to_fine_cells[level].shape[1] * level_ratio + key = (("inject", ncandidate) + (Vf.block_size,) + entity_dofs_key(Vc.finat_element.complex.get_topology()) + entity_dofs_key(Vf.finat_element.complex.get_topology()) @@ -364,7 +365,6 @@ def inject_kernel(Vf, Vc): try: return cache[key] except KeyError: - ncandidate = hierarchy.coarse_to_fine_cells[level].shape[1] * level_ratio return cache.setdefault(key, (dg_injection_kernel(Vf, Vc, ncandidate), True)) else: expression = ufl.Coefficient(Vf) @@ -403,6 +403,7 @@ def set_coordinates(self, domain): self._coefficient(f, "macro_coords") def _coefficient(self, coefficient, name): + """Register a coefficient as a macro-cell kernel argument and return its GEM expression.""" element = create_element(coefficient.ufl_element()) shape = self.shape + element.index_shape size = numpy.prod(shape, dtype=int) diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index fecbae149e..8c1a6c237a 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -1,4 +1,5 @@ import numpy as np +import warnings from fractions import Fraction from collections import defaultdict from collections.abc import Sequence @@ -7,40 +8,60 @@ import petsctools import firedrake -import firedrake.cython.dmcommon as dmcommon from functools import cached_property from firedrake import utils from firedrake.cython import mgimpl as impl +import firedrake.cython.dmcommon as dmcommon from .utils import set_level __all__ = ("HierarchyBase", "MeshHierarchy", "ExtrudedMeshHierarchy", "NonNestedHierarchy", "SemiCoarsenedExtrudedHierarchy", "SubmeshHierarchy") +def make_unoverlapped_dm(dm): + """Effectively invert dm.distributeOverlap(). + + The resulting plex has the identical data structure as the one before + distributeOverlap(). This is algorithmically guaranteed. + """ + tdim = dm.getDimension() + dm = dmcommon.submesh_create(dm, tdim, "depth", tdim, ignore_label_halo=True) + dm.removeLabel("pyop2_core") + dm.removeLabel("pyop2_owned") + dm.removeLabel("pyop2_ghost") + return dm + + class HierarchyBase(object): """Create an encapsulation of an hierarchy of meshes. - :arg meshes: list of meshes (coarse to fine) - :arg coarse_to_fine_cells: list of numpy arrays for each level - pair, mapping each coarse cell into fine cells it intersects. - :arg fine_to_coarse_cells: list of numpy arrays for each level - pair, mapping each fine cell into coarse cells it intersects. - :arg refinements_per_level: number of mesh refinements each - multigrid level should "see". - :arg nested: Is this mesh hierarchy nested? - - .. note:: + Parameters + ---------- + meshes : + List of meshes (coarse to fine). + coarse_to_fine_cells : + List of numpy arrays for each level pair, mapping each coarse cell + into fine cells it intersects. + fine_to_coarse_cells : + List of numpy arrays for each level pair, mapping each fine cell into + coarse cells it intersects. + refinements_per_level : + Number of mesh refinements each multigrid level should "see". + nested : + Is this mesh hierarchy nested? + + Notes + ----- + Most of the time, you do not need to create this object yourself, instead + using `MeshHierarchy`, `ExtrudedMeshHierarchy`, or `NonNestedHierarchy`. - Most of the time, you do not need to create this object - yourself, instead using :func:`MeshHierarchy`, - :func:`ExtrudedMeshHierarchy`, or :func:`NonNestedHierarchy`. """ def __init__(self, meshes, coarse_to_fine_cells, fine_to_coarse_cells, refinements_per_level=1, nested=False): petsctools.cite("Mitchell2016") - self._meshes = tuple(meshes) - self.meshes = tuple(meshes[::refinements_per_level]) + self._meshes = list(meshes) + self.meshes = self._meshes[::refinements_per_level] self.coarse_to_fine_cells = coarse_to_fine_cells self.fine_to_coarse_cells = fine_to_coarse_cells self.refinements_per_level = refinements_per_level @@ -73,13 +94,99 @@ def __getitem__(self, idx): :arg idx: The :func:`~.Mesh` to return""" return self.meshes[idx] - -def MeshHierarchy(mesh, refinement_levels, + def add_mesh(self, mesh, coarse_to_fine_cells=None, fine_to_coarse_cells=None): + """Add a mesh on top of the finest level of the hierarchy. + + Only supported for hierarchies with ``refinements_per_level == 1``. + + Parameters + ---------- + mesh : + The mesh to add, usually obtained by calling + :meth:`~firedrake.mesh.MeshGeometry.refine_marked_elements` on the + current finest mesh. + coarse_to_fine_cells : + Map from the cells of the current finest mesh to the cells of + ``mesh``. Defaults to the map ``mesh`` recorded when it was + adaptively refined. + fine_to_coarse_cells : + Map from the cells of ``mesh`` to the cells of the current finest + mesh. Defaults the same way as ``coarse_to_fine_cells``. + + Returns + ------- + MeshGeometry + The mesh that was added. + + """ + if self.refinements_per_level != 1: + raise NotImplementedError("Cannot add a mesh to a hierarchy with " + "refinements_per_level > 1") + if coarse_to_fine_cells is None or fine_to_coarse_cells is None: + if mesh.adaptive_parent is self[-1]: + coarse_to_fine_cells, fine_to_coarse_cells = mesh.adaptive_cell_maps + elif self.nested: + raise ValueError("Expecting a mesh adaptively refined from the finest " + "level of this hierarchy, or explicit cell maps") + + level = len(self.meshes) + self._meshes.append(mesh) + self.meshes.append(mesh) + set_level(mesh, self, level) + mesh.topology_dm.setRefineLevel(level) + self.coarse_to_fine_cells[Fraction(level - 1, 1)] = coarse_to_fine_cells + self.fine_to_coarse_cells[Fraction(level, 1)] = fine_to_coarse_cells + return mesh + + def adapt(self, eta, theta: float): + """Add a new mesh to the hierarchy by locally refining the finest mesh + with a simplified variant of Dorfler marking. + + Parameters + ---------- + eta : + A DG0 :class:`~firedrake.function.Function` with the local error estimator. + theta : + The threshold for marking as a fraction of the maximum error. + + Returns + ------- + MeshGeometry + The mesh that was added. + + Note + ---- + Dorfler marking involves sorting all of the elements by decreasing + error estimator and taking the minimal set that exceeds some fixed + fraction of the total error. What this code implements is the simpler + variant that doesn't have a proof of convergence (as far as I know) + but works as well in practice. + + """ + if not isinstance(eta, (firedrake.Function, firedrake.Cofunction)): + raise TypeError(f"eta must be a Function or Cofunction, not a {type(eta).__name__}") + M = eta.function_space() + if M.finat_element.space_dimension() != 1: + raise ValueError("eta must be a Function or Cofunction in DG0") + mesh = self[-1] + if M.mesh() is not mesh: + raise ValueError("eta must be defined on the finest mesh of the hierarchy") + + # Take the maximum over all processes + with eta.dat.vec_ro as evec: + _, eta_max = evec.max() + + markers = firedrake.Function(M) + markers.dat.data_wo[eta.dat.data_ro > theta * eta_max] = 1 + return self.add_mesh(mesh.refine_marked_elements(markers)) + + +def MeshHierarchy(mesh, refinement_levels=0, refinements_per_level=1, netgen_flags=False, reorder=None, distribution_parameters=None, callbacks=None, - mesh_builder=firedrake.Mesh): + mesh_builder=firedrake.Mesh, nested=True): """Build a hierarchy of meshes by uniformly refining a coarse mesh. Parameters @@ -87,9 +194,12 @@ def MeshHierarchy(mesh, refinement_levels, mesh : MeshGeometry the coarse mesh to refine refinement_levels : int - the number of levels of refinement + the number of levels of uniform refinement. This may be dynamically + increased by :meth:`HierarchyBase.adapt` or + :meth:`HierarchyBase.add_mesh`. refinements_per_level : int the number of refinements for each level in the hierarchy. + Adaptive refinement only supports one refinement per level. netgen_flags : bool, dict either a bool or a dictionary containing options for Netgen. If not False the hierachy is constructed using ngsPETSc, if @@ -109,10 +219,16 @@ def MeshHierarchy(mesh, refinement_levels, callback receives the refined DM (and the current level). mesh_builder Function to turn a DM into a ``Mesh``. Used by pyadjoint. + nested : bool + Are the meshes added to this hierarchy required to be nested? If + `False`, :meth:`HierarchyBase.add_mesh` accepts a mesh that was not + adaptively refined from the finest level. + Returns ------- - A :py:class:`HierarchyBase` object representing the - mesh hierarchy. + HierarchyBase + The mesh hierarchy. + """ if (isinstance(netgen_flags, bool) and netgen_flags) or isinstance(netgen_flags, dict): @@ -124,60 +240,62 @@ def MeshHierarchy(mesh, refinement_levels, else: raise RuntimeError("Cannot create a NetgenHierarchy from a mesh that has not been generated by\ Netgen.") - # Effectively "invert" addOverlap(). - # -- The resulting plex is to have the identical data structure as the one before addOverlap(). - # This is algorithmically guaranteed. - tdim = mesh.topology_dm.getDimension() - cdm = dmcommon.submesh_create(mesh.topology_dm, tdim, "depth", tdim, True) - cdm.removeLabel("pyop2_core") - cdm.removeLabel("pyop2_owned") - cdm.removeLabel("pyop2_ghost") - cdm.setRefinementUniform(True) - dms = [cdm] if callbacks is not None: before, after = callbacks else: before = after = lambda dm, i: None + + # Refine an unoverlapped plex at each level. Keeping every dm here + # unoverlapped means overlap only ever needs to be added once, by + # mesh_builder below. + cdm = mesh.topology_dm + if refinement_levels > 0: + cdm = make_unoverlapped_dm(cdm) + cdm.setRefinementUniform(True) + dms = [cdm] for i in range(refinement_levels*refinements_per_level): if i % refinements_per_level == 0: before(cdm, i) rdm = cdm.refine() if i % refinements_per_level == 0: after(rdm, i) - dms.append(rdm) - cdm = rdm # Fix up coords if refining embedded circle or sphere if hasattr(mesh, '_radius'): # FIXME, really we need some CAD-like representation # of the boundary we're trying to conform to. This # doesn't DTRT really for cubed sphere meshes (the # refined meshes are no longer gnonomic). - coords = cdm.getCoordinatesLocal().array.reshape(-1, mesh.geometric_dimension) + coords = rdm.getCoordinatesLocal().array.reshape(-1, mesh.geometric_dimension) scale = mesh._radius / np.linalg.norm(coords, axis=1).reshape(-1, 1) coords *= scale - lgmaps_without_overlap = [impl.create_lgmap(dm) for dm in dms] + + dms.append(rdm) + cdm = rdm + + # Build a mesh for each level, adding overlap here. parameters = {} if distribution_parameters is not None: parameters.update(distribution_parameters) else: parameters.update(mesh._distribution_parameters) parameters["partition"] = False - meshes = [mesh] + [ - mesh_builder( - dm, + + meshes = [mesh] + for rdm in dms[1:]: + fmesh = mesh_builder( + rdm, dim=mesh.geometric_dimension, distribution_parameters=parameters, reorder=reorder, comm=mesh.comm, ) - for dm in dms[1:] - ] - lgmaps_with_overlap = [] - for i, m in enumerate(meshes): - lgmaps_with_overlap.append(impl.create_lgmap(m.topology_dm)) - m.topology_dm.setRefineLevel(i) + meshes.append(fmesh) + + # Build local-to-global maps and coarse/fine cell maps between + # consecutive levels. lgmaps = [ - (no, o) for no, o in zip(lgmaps_without_overlap, lgmaps_with_overlap) + (impl.create_lgmap(dm), impl.create_lgmap(m.topology_dm)) + for dm, m in zip(dms, meshes) ] coarse_to_fine_cells = [] fine_to_coarse_cells = [None] @@ -187,12 +305,16 @@ def MeshHierarchy(mesh, refinement_levels, coarse_to_fine_cells.append(c2f) fine_to_coarse_cells.append(f2c) + for i, m in enumerate(meshes): + # Firedrake counts multigrid levels, PETSc counts refinements + m.topology_dm.setRefineLevel(i) + coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), c2f) for i, c2f in enumerate(coarse_to_fine_cells)) fine_to_coarse_cells = dict((Fraction(i, refinements_per_level), f2c) for i, f2c in enumerate(fine_to_coarse_cells)) return HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, - refinements_per_level, nested=True) + refinements_per_level, nested=nested) def ExtrudedMeshHierarchy(base_hierarchy, height, base_layer=-1, refinement_ratio=2, layers=None, @@ -308,6 +430,15 @@ def SemiCoarsenedExtrudedHierarchy(base_mesh, height, nref=1, base_layer=-1, ref nested=True) +def AdaptiveMeshHierarchy(*args, **kwargs): + """Deprecated alias for `MeshHierarchy`.""" + warnings.warn( + "The ``AdaptiveMeshHierarchy`` constructor is deprecated and will be removed in a future " + "release. Please use the ``MeshHierarchy`` constructor instead.", FutureWarning + ) + return MeshHierarchy(*args, **kwargs) + + def NonNestedHierarchy(*meshes): return HierarchyBase(meshes, [None for _ in meshes], [None for _ in meshes], nested=False) diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index f6088528ed..a957e50c7d 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -278,7 +278,9 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): comm = mesh.comm for l in range(1, levs+1): rdm, ngmesh = refinementTypes[refType][0](base_ngmesh, cdm) - cdm = rdm + # `fd.Mesh` mutates `rdm` in place (e.g. adding overlap), so clone + # it first to keep an unoverlapped dm for the next refinement. + cdm = rdm.clone() if optMoves: # Optimises the mesh, for example smoothing if tdim == 2: @@ -349,5 +351,6 @@ def reconstruct_mesh(mesh, *args, **kwargs): tmesh._did_reordering = mesh._did_reordering tmesh.netgen_mesh = mesh.netgen_mesh tmesh.netgen_flags = mesh.netgen_flags + tmesh.sfBC = mesh.sfBC tmesh.sfBC_orig = mesh.sfBC_orig return tmesh diff --git a/firedrake/mg/ufl_utils.py b/firedrake/mg/ufl_utils.py index e5f1972cdc..f7eb888fcf 100644 --- a/firedrake/mg/ufl_utils.py +++ b/firedrake/mg/ufl_utils.py @@ -71,6 +71,7 @@ def facet_normal(self, o): @singledispatch def _reconstruct(expr, self, coefficient_mapping=None): + """Fallback case: leave an expression with no registered handler unchanged.""" return expr @@ -559,6 +560,14 @@ def mult(self, mat, x, y): with self.cfn.dat.vec_ro as v: v.copy(y) + def multTranspose(self, mat, x, y): + # PETSc's MatRestrict() cannot distinguish an injection matrix from + # an interpolation matrix when the coarse and fine spaces happen to + # have equal size (e.g. a no-op adaptive refinement level), and may + # call MatMultTranspose() instead of MatMult(). Injection is only + # ever used in one direction (fine to coarse), so both must agree. + self.mult(mat, x, y) + def create_interpolation(dmc, dmf): cctx = get_appctx(dmc) diff --git a/firedrake/mg/utils.py b/firedrake/mg/utils.py index d2c37ed7aa..ffc3c1e8ad 100644 --- a/firedrake/mg/utils.py +++ b/firedrake/mg/utils.py @@ -73,6 +73,27 @@ def coarse_node_to_fine_node_map(Vc, Vf): coarse_to_fine = hierarchy.coarse_to_fine_cells[levelc] coarse_to_fine_nodes = impl.coarse_to_fine_nodes(Vc, Vf, coarse_to_fine) + # Under adaptive refinement, coarse cells have varying numbers of + # fine descendants, so coarse_to_fine (and hence coarse_to_fine_nodes) + # is right-padded with -1 up to the busiest coarse cell's count. + # op2.Map cannot hold negative indices, and every *owned* coarse + # node needs at least one real candidate to inject from; but padding + # slots on rows that do have candidates can safely be filled with a + # duplicate of one of that row's real entries; the injection kernel + # below only ever reads (op2.READ) through this map and picks the + # candidate matching the coarse node's physical location, so a + # repeated valid entry is just redundantly (harmlessly) considered. + valid = coarse_to_fine_nodes >= 0 + if not valid.all(): + nonempty = valid.any(axis=1) + if not nonempty[:Vc.node_set.size].all(): + raise RuntimeError("Adaptive coarse-to-fine map has empty node candidates") + replacement = numpy.zeros(coarse_to_fine_nodes.shape[0], + dtype=coarse_to_fine_nodes.dtype) + rows = numpy.nonzero(nonempty)[0] + replacement[rows] = coarse_to_fine_nodes[rows, valid[rows].argmax(axis=1)] + coarse_to_fine_nodes = numpy.where(valid, coarse_to_fine_nodes, + replacement[:, None]) return cache.setdefault(key, op2.Map(Vc.node_set, Vf.node_set, coarse_to_fine_nodes.shape[1], values=coarse_to_fine_nodes)) @@ -110,13 +131,20 @@ def coarse_cell_to_fine_node_map(Vc, Vf): coarse_to_fine = hierarchy.coarse_to_fine_cells[levelc] _, ncell = coarse_to_fine.shape iterset = Vc.mesh().cell_set - arity = Vf.finat_element.space_dimension() * ncell + fine_per_cell = Vf.finat_element.space_dimension() + arity = fine_per_cell * ncell coarse_to_fine_nodes = numpy.full((iterset.total_size, arity*level_ratio), -1, dtype=IntType) - values = Vf.cell_node_map().values[coarse_to_fine, :].reshape(iterset.size, arity) + values = numpy.full((iterset.size, ncell, fine_per_cell), -1, dtype=IntType) + owned_coarse_to_fine = coarse_to_fine[:iterset.size, :] + valid = owned_coarse_to_fine >= 0 + values[valid, :] = Vf.cell_node_map().values[owned_coarse_to_fine[valid], :] + values = values.reshape(iterset.size, arity) if Vc.extruded: off = numpy.tile(Vf.offset, ncell) - coarse_to_fine_nodes[:Vc.mesh().cell_set.size, :] = numpy.hstack([values + off*i for i in range(level_ratio)]) + coarse_to_fine_nodes[:Vc.mesh().cell_set.size, :] = numpy.hstack([ + numpy.where(values >= 0, values + off*i, -1) for i in range(level_ratio) + ]) else: coarse_to_fine_nodes[:Vc.mesh().cell_set.size, :] = values offset = Vf.offset diff --git a/firedrake/netgen.py b/firedrake/netgen.py index e44af74140..c0c0369779 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -15,7 +15,7 @@ try: import netgen.meshing as ngm from netgen.meshing import MeshingParameters - from ngsPETSc import MeshMapping + from ngsPETSc import MeshMapping, createNetgenMesh except ImportError: pass @@ -129,6 +129,42 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray): return permutation +def _transfer_high_order_coordinates(coarse_mesh, fine_mesh, order): + """Transfer high-order coordinates from a Netgen geometry to a refined mesh. + + ``fine_mesh`` is a straight-edged (order 1) refinement of ``coarse_mesh``. + This rebuilds its Netgen mesh from ``coarse_mesh``'s geometry and curves + it to the requested ``order``, so that the curved fine mesh follows the + same underlying CAD geometry as the coarse one, rather than just + interpolating the coarse mesh's straight-edged coordinates. + + Parameters + ---------- + coarse_mesh : MeshGeometry + The coarse mesh, carrying the Netgen geometry to curve against. + fine_mesh : MeshGeometry + A straight-edged refinement of ``coarse_mesh``. Its Netgen attributes + are set here, as they are required to curve it. + order : int + The polynomial order of the curved coordinate field. + + Returns + ------- + MeshGeometry + A mesh sharing ``fine_mesh``'s topology, with coordinates curved to + ``order`` against ``coarse_mesh``'s geometry. + + """ + fine_mesh.netgen_mesh = createNetgenMesh(fine_mesh.topology_dm, coarse_mesh.netgen_mesh) + fine_mesh.netgen_flags = getattr(coarse_mesh, "netgen_flags", {}) + cg_field = not coarse_mesh.coordinates.function_space().finat_element.is_dg() + curved_coordinates = fine_mesh.curve_field(order=order, cg_field=cg_field) + curved_mesh = firedrake.Mesh(curved_coordinates, name=fine_mesh.name) + curved_mesh.netgen_mesh = fine_mesh.netgen_mesh + curved_mesh.netgen_flags = fine_mesh.netgen_flags + return curved_mesh + + def splitToQuads(plex, dim, comm): """Split a Netgen mesh into quads using a PETSc transform.""" # TODO: Improve support quad meshing. diff --git a/pyproject.toml b/pyproject.toml index 063133c4e1..1aa437ccd1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,7 +33,8 @@ dependencies = [ "packaging", # TODO RELEASE # "petsc4py==3.25.0", - "petsctools>=2026.0", + # TODO RELEASE + "petsctools @ git+https://github.com/firedrakeproject/petsctools.git@main", "pkgconfig", "progress", "pyadjoint-ad>=2026.4.0", diff --git a/tests/firedrake/multigrid/test_adaptive_multigrid.py b/tests/firedrake/multigrid/test_adaptive_multigrid.py index 1235c95bd5..e333c3be74 100644 --- a/tests/firedrake/multigrid/test_adaptive_multigrid.py +++ b/tests/firedrake/multigrid/test_adaptive_multigrid.py @@ -1,106 +1,308 @@ -""" -Tests for AdaptiveMeshHierarchy -and AdaptiveTransferManager -""" - import pytest import numpy as np +from mpi4py import MPI from firedrake import * -@pytest.fixture(params=[2, 3]) -def amh(request): - """ - Generate AdaptiveMeshHierarchies - """ - from netgen.occ import WorkPlane, OCCGeometry, Box, Pnt - dim = request.param - if dim == 2: +def corner_adaptive_hierarchy(base, nlevels): + mh = MeshHierarchy(base) + for l in range(nlevels): + mesh = mh[-1] + x = SpatialCoordinate(mesh) + M = FunctionSpace(mesh, "DG", 0) + m = Function(M, name="marker") + m.interpolate(conditional(sum(x) < 2**(-l+1), 1, 0)) + mh.add_mesh(mesh.refine_marked_elements(m)) + return mh + + +def _linear_expr(mesh): + """A linear expression in the mesh's spatial coordinates, generalizing + ``x + 2*y`` to any dimension (``x + 2*y + 3*z`` in 3D, etc.).""" + x = SpatialCoordinate(mesh) + weights = Constant(list(range(1, mesh.geometric_dimension + 1))) + return dot(weights, x) + + +@pytest.fixture(params=[ + "firedrake-square", + "firedrake-cube", + pytest.param("netgen-square", marks=pytest.mark.skipnetgen), + pytest.param("netgen-cube", marks=pytest.mark.skipnetgen), +]) +def coarse_mesh(request): + dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)} + mesher = request.param + if mesher == "firedrake-square": + return UnitSquareMesh(1, 1, distribution_parameters=dparams) + elif mesher == "firedrake-cube": + return UnitCubeMesh(1, 1, 1, distribution_parameters=dparams) + elif mesher == "netgen-square": + from netgen.occ import WorkPlane, OCCGeometry wp = WorkPlane() wp.Rectangle(1, 1) face = wp.Face() geo = OCCGeometry(face, dim=2) - maxh = 0.5 - else: + ngmesh = geo.GenerateMesh(maxh=0.5) + return Mesh(ngmesh, distribution_parameters=dparams) + elif mesher == "netgen-cube": + from netgen.occ import Box, OCCGeometry, Pnt cube = Box(Pnt(0, 0, 0), Pnt(1, 1, 1)) geo = OCCGeometry(cube, dim=3) - maxh = 0.5 + ngmesh = geo.GenerateMesh(maxh=0.5) + return Mesh(ngmesh, distribution_parameters=dparams) + else: + raise NotImplementedError(f"Unrecognized mesher {mesher}") - ngmesh = geo.GenerateMesh(maxh=maxh) - dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)} - base = Mesh(ngmesh, distribution_parameters=dparams) - amh_test = AdaptiveMeshHierarchy(base) +@pytest.fixture +def mh(coarse_mesh): + return corner_adaptive_hierarchy(coarse_mesh, nlevels=2) - rg = RandomGenerator(PCG64(seed=0)) - for l in range(2): - mesh = amh_test[-1] - DG = FunctionSpace(mesh, "DG", 0) - should_refine = rg.uniform(DG, 0, 1).dat.global_data - ngmesh = mesh.netgen_mesh - if dim == 2: - els = ngmesh.Elements2D() - else: - els = ngmesh.Elements3D() - for i, el in enumerate(els): - el.refine = 1 if should_refine[i] < 0.5 else 0 +def test_refine_marked_elements_populates_cell_maps(coarse_mesh): + mesh = coarse_mesh + mh = MeshHierarchy(mesh) - ngmesh.Refine(adaptive=True) - mesh = Mesh(ngmesh, distribution_parameters=dparams) - amh_test.add_mesh(mesh) - return amh_test + M = FunctionSpace(mesh, "DG", 0) + markers = Function(M) + markers.dat.data_wo[0] = 1 + refined_mesh = mesh.refine_marked_elements(markers) + mh.add_mesh(refined_mesh) -@pytest.fixture -def mh_uniform(): - """ - Generate MeshHierarchy for reference - """ - from netgen.occ import WorkPlane, OCCGeometry - wp = WorkPlane() - wp.Rectangle(2, 2) - face = wp.Face() - geo = OCCGeometry(face, dim=2) - maxh = 0.5 - ngmesh = geo.GenerateMesh(maxh=maxh) + coarse_to_fine = mh.coarse_to_fine_cells[0] + fine_to_coarse = mh.fine_to_coarse_cells[1] - dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)} - base1 = Mesh(ngmesh, distribution_parameters=dparams) - mh = MeshHierarchy(base1, 2) + assert coarse_to_fine.shape[0] == mesh.cell_set.size + assert fine_to_coarse.shape == (refined_mesh.cell_set.size, 1) + assert (fine_to_coarse >= -1).all() + assert (fine_to_coarse >= 0).any() + assert (coarse_to_fine >= 0).any() + for coarse_cell, fine_cells in enumerate(coarse_to_fine): + fine_cells = fine_cells[(fine_cells >= 0) & (fine_cells < fine_to_coarse.shape[0])] + if fine_cells.size: + assert (fine_to_coarse[fine_cells, 0] == coarse_cell).all() - base2 = Mesh(ngmesh, distribution_parameters=dparams) - amh = AdaptiveMeshHierarchy(base2) - for _ in range(2): - mesh = amh[-1] - ngmesh = mesh.netgen_mesh - ngmesh.Refine() - mesh = Mesh(ngmesh, distribution_parameters=dparams) - amh.add_mesh(mesh) - return amh, mh +def test_refine_marked_elements_is_local(): + # Regression test: dm_plex_transform_type=refine_sbr must actually reach + # PETSc's options database. If it doesn't, dm.adaptLabel() silently + # falls back to unconditional uniform refinement of every cell, + # regardless of marking. + nx = 8 + mesh = UnitSquareMesh(nx, nx) + ncoarse = mesh.cell_set.size -@pytest.fixture -def atm(): - """atm used in tests""" - return AdaptiveTransferManager() + M = FunctionSpace(mesh, "DG", 0) + markers = Function(M) + markers.dat.data_wo[0] = 1 + refined_mesh = mesh.refine_marked_elements(markers) + coarse_to_fine, _ = refined_mesh.adaptive_cell_maps -@pytest.fixture -def tm(): - """tm used for restrict consistency""" - return TransferManager() + n_children = (coarse_to_fine >= 0).sum(axis=1) + unmarked = np.ones(ncoarse, dtype=bool) + unmarked[0] = False + + # refine_sbr's conforming closure may also split a handful of + # neighbouring cells (to avoid hanging nodes), but with a marked region + # this small, most of the mesh must be left untouched. + assert (n_children[unmarked] == 1).sum() >= 0.5 * unmarked.sum() @pytest.mark.parallel([1, 2]) +def test_refine_marked_elements_repeats(coarse_mesh): + """A marker value of n refines the marked cells n times, and the cell maps + reach all the way from the original mesh to the n-times-refined one.""" + mesh = coarse_mesh + ncells = {} + max_children = {} + for n in (1, 2): + M = FunctionSpace(mesh, "DG", 0) + markers = Function(M) + markers.dat.data_wo[:1] = n + + refined_mesh = mesh.refine_marked_elements(markers) + coarse_to_fine, fine_to_coarse = refined_mesh.adaptive_cell_maps + + assert coarse_to_fine.shape[0] == mesh.cell_set.size + assert fine_to_coarse.shape == (refined_mesh.cell_set.size, 1) + for coarse_cell, fine_cells in enumerate(coarse_to_fine): + fine_cells = fine_cells[(fine_cells >= 0) & (fine_cells < fine_to_coarse.shape[0])] + assert (fine_to_coarse[fine_cells, 0] == coarse_cell).all() + assert np.allclose(assemble(1*dx(refined_mesh)), assemble(1*dx(mesh))) + + ncells[n] = mesh.comm.allreduce(refined_mesh.cell_set.size) + max_children[n] = mesh.comm.allreduce( + (coarse_to_fine >= 0).sum(axis=1).max(initial=0), op=MPI.MAX) + + # Each extra round bisects the marked cells again, so they gain both cells + # overall and descendants of the cell they came from. + assert ncells[2] > ncells[1] + assert max_children[2] > max_children[1] + + +def test_add_mesh_rejects_unrelated_mesh(): + """Cell maps are only meaningful relative to the mesh they were built + against, so a mesh refined from anything but the finest level is refused + rather than silently recorded with somebody else's maps.""" + mh = MeshHierarchy(UnitSquareMesh(2, 2)) + + other = UnitSquareMesh(4, 4) + assert other.adaptive_parent is None + with pytest.raises(ValueError): + mh.add_mesh(other) + + markers = Function(FunctionSpace(other, "DG", 0)) + markers.dat.data_wo[:1] = 1 + foreign = other.refine_marked_elements(markers) + assert foreign.adaptive_parent is other + with pytest.raises(ValueError): + mh.add_mesh(foreign) + + +@pytest.mark.parallel([1, 2, 4]) +def test_adapt_basic(): + nx = 1 + base = UnitCubeMesh(nx, nx, nx) + + mh = corner_adaptive_hierarchy(base, nlevels=6) + + mesh = mh[-1] + assert np.allclose(assemble(1*dx(mesh)), assemble(1*dx(base))) + + +def test_CG1_native_transfers_use_adaptive_cell_maps(coarse_mesh): + mesh = coarse_mesh + mh = MeshHierarchy(mesh) + + M = FunctionSpace(mesh, "DG", 0) + markers = Function(M) + markers.dat.data_wo[0] = 1 + refined_mesh = mesh.refine_marked_elements(markers) + mh.add_mesh(refined_mesh) + + assert (mh.coarse_to_fine_cells[0] < 0).any() + + V_coarse = FunctionSpace(mesh, "CG", 1) + V_fine = FunctionSpace(refined_mesh, "CG", 1) + expr_coarse = _linear_expr(mesh) + expr_fine = _linear_expr(refined_mesh) + + u_coarse = Function(V_coarse).interpolate(expr_coarse) + u_fine = Function(V_fine) + prolong(u_coarse, u_fine) + assert errornorm(expr_fine, u_fine) <= 1e-12 + + u_fine_exact = Function(V_fine).interpolate(expr_fine) + u_coarse_injected = Function(V_coarse) + inject(u_fine_exact, u_coarse_injected) + assert errornorm(expr_coarse, u_coarse_injected) <= 1e-12 + + r_fine = assemble(conj(TestFunction(V_fine)) * dx) + r_coarse = Cofunction(V_coarse.dual()) + restrict(r_fine, r_coarse) + assert np.allclose( + assemble(action(r_coarse, u_coarse)), + assemble(action(r_fine, u_fine)), + rtol=1e-12, + atol=1e-12, + ) + + +def _assert_adapt_after_uniform_refinement(mh): + """Adaptively refine the finest level of the uniformly-refined hierarchy + ``mh`` by marking a single cell, and check that the cell maps of the level + this adds are sane. Shared by the ``test_adapt_after_uniform_*refinement`` + tests, which only differ in how ``mh`` itself was built. + """ + mesh = mh[-1] + level = len(mh) + + M = FunctionSpace(mesh, "DG", 0) + markers = Function(M) + # Slice rather than index: with more ranks than coarse cells, some ranks + # legitimately own zero local cells, and an unconditional [0] = 1 would + # raise IndexError there, desynchronizing the collectives inside + # refine_marked_elements and hanging the surviving ranks. + markers.dat.data_wo[:1] = 1 + + refined_mesh = mh.add_mesh(mesh.refine_marked_elements(markers)) + assert len(mh) == level + 1 + assert mh[-1] is refined_mesh + + coarse_to_fine = mh.coarse_to_fine_cells[level - 1] + fine_to_coarse = mh.fine_to_coarse_cells[level] + + assert coarse_to_fine.shape[0] == mesh.cell_set.size + assert fine_to_coarse.shape == (refined_mesh.cell_set.size, 1) + # A rank may legitimately own zero local cells (e.g. more ranks than + # coarse cells), leaving these arrays empty on that rank alone, so the + # "some entry is valid" check must be collective, not per-rank. + assert mesh.comm.allreduce((fine_to_coarse >= 0).any(), op=MPI.LOR) + assert mesh.comm.allreduce((coarse_to_fine >= 0).any(), op=MPI.LOR) + + @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 2, 4]) +def test_adapt_after_uniform_netgen_refinement(): + from netgen.geom2d import SplineGeometry + + geo = SplineGeometry() + geo.AddRectangle((0, 0), (1, 1), bc="boundary") + netgen_mesh = geo.GenerateMesh(maxh=0.5) + netgen_mesh.Refine() + mesh = Mesh(netgen_mesh) + _assert_adapt_after_uniform_refinement(MeshHierarchy(mesh)) + + +@pytest.mark.skipnetgen +@pytest.mark.parallel([1, 2]) +@pytest.mark.parametrize("degree", [1, 2]) +def test_adapt_preserves_mesh_metadata(degree): + """Adaptive refinement carries the Netgen geometry and flags, and the mesh + construction parameters, over to the refined mesh, so that the refined + mesh can itself be refined again.""" + from netgen.geom2d import CSG2d, Circle + geo = CSG2d() + geo.Add(Circle(center=(0, 0), radius=1.0, bc="circle")) + ngmesh = geo.GenerateMesh(maxh=0.75) + dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 2)} + mesh = Mesh(ngmesh, netgen_flags={"degree": degree}, + distribution_parameters=dparams) + + markers = Function(FunctionSpace(mesh, "DG", 0)).assign(1) + refined = mesh.refine_marked_elements(markers) + + assert refined.netgen_mesh is not None + assert refined.netgen_flags == mesh.netgen_flags + assert refined._distribution_parameters == mesh._distribution_parameters + assert refined.tolerance == mesh.tolerance + assert refined.coordinates.function_space().ufl_element().degree() == degree + + markers = Function(FunctionSpace(refined, "DG", 0)).assign(1) + twice_refined = refined.refine_marked_elements(markers) + assert twice_refined.netgen_flags == mesh.netgen_flags + assert twice_refined._distribution_parameters == mesh._distribution_parameters + assert twice_refined.coordinates.function_space().ufl_element().degree() == degree + + +@pytest.mark.parallel([1, 2, 4]) +@pytest.mark.parametrize("refine", [1, 2]) +def test_adapt_after_uniform_refinement(coarse_mesh, refine): + """A hierarchy built by uniform refinement can be adaptively refined.""" + netgen_flags = {} if hasattr(coarse_mesh, "netgen_mesh") else None + mh = MeshHierarchy(coarse_mesh, refine, netgen_flags=netgen_flags) + _assert_adapt_after_uniform_refinement(mh) + + +@pytest.mark.parallel([1, 2, 4]) @pytest.mark.parametrize("operator", ["prolong", "inject"]) -def test_DG0(amh, atm, operator): # pylint: disable=W0621 - """ - Prolongation & Injection test for DG0 - """ - V_coarse = FunctionSpace(amh[0], "DG", 0) - V_fine = FunctionSpace(amh[-1], "DG", 0) +def test_DG0(mh, operator): + """Prolongation & Injection test for DG0""" + V_coarse = FunctionSpace(mh[0], "DG", 0) + V_fine = FunctionSpace(mh[-1], "DG", 0) u_coarse = Function(V_coarse) u_fine = Function(V_fine) xc, *_ = SpatialCoordinate(V_coarse.mesh()) @@ -112,25 +314,22 @@ def test_DG0(amh, atm, operator): # pylint: disable=W0621 u_coarse.interpolate(stepc) assert errornorm(stepc, u_coarse) <= 1e-12 - atm.prolong(u_coarse, u_fine) + prolong(u_coarse, u_fine) assert errornorm(stepf, u_fine) <= 1e-12 if operator == "inject": u_fine.interpolate(stepf) assert errornorm(stepf, u_fine) <= 1e-12 - atm.inject(u_fine, u_coarse) + inject(u_fine, u_coarse) assert errornorm(stepc, u_coarse) <= 1e-12 -@pytest.mark.parallel([1, 2]) -@pytest.mark.skipnetgen +@pytest.mark.parallel([1, 2, 4]) @pytest.mark.parametrize("operator", ["prolong", "inject"]) -def test_CG1(amh, atm, operator): # pylint: disable=W0621 - """ - Prolongation & Injection test for CG1 - """ - V_coarse = FunctionSpace(amh[0], "CG", 1) - V_fine = FunctionSpace(amh[-1], "CG", 1) +def test_CG1(mh, operator): + """Prolongation & Injection test for CG1""" + V_coarse = FunctionSpace(mh[0], "CG", 1) + V_fine = FunctionSpace(mh[-1], "CG", 1) u_coarse = Function(V_coarse) u_fine = Function(V_fine) xc, *_ = SpatialCoordinate(V_coarse.mesh()) @@ -140,80 +339,31 @@ def test_CG1(amh, atm, operator): # pylint: disable=W0621 u_coarse.interpolate(xc) assert errornorm(xc, u_coarse) <= 1e-12 - atm.prolong(u_coarse, u_fine) + prolong(u_coarse, u_fine) assert errornorm(xf, u_fine) <= 1e-12 if operator == "inject": u_fine.interpolate(xf) assert errornorm(xf, u_fine) <= 1e-12 - atm.inject(u_fine, u_coarse) + inject(u_fine, u_coarse) assert errornorm(xc, u_coarse) <= 1e-12 -@pytest.mark.parallel([1, 2]) -@pytest.mark.skipnetgen -def test_restrict_consistency(mh_uniform, atm, tm): # pylint: disable=W0621 - """ - Test restriction consistency of amh with uniform refinement vs mh - """ - amh_unif, mh = mh_uniform - - V_coarse = FunctionSpace(amh_unif[0], "DG", 0) - V_fine = FunctionSpace(amh_unif[-1], "DG", 0) - u_coarse = Function(V_coarse) - u_fine = Function(V_fine) - xc, _ = SpatialCoordinate(V_coarse.mesh()) - - u_coarse.interpolate(xc) - atm.prolong(u_coarse, u_fine) - - rf = assemble(conj(TestFunction(V_fine)) * dx) - rc = Cofunction(V_coarse.dual()) - atm.restrict(rf, rc) - - # compare with mesh_hierarchy - xcoarse, _ = SpatialCoordinate(mh[0]) - Vcoarse = FunctionSpace(mh[0], "DG", 0) - Vfine = FunctionSpace(mh[-1], "DG", 0) - - mhuc = Function(Vcoarse) - mhuc.interpolate(xcoarse) - mhuf = Function(Vfine) - tm.prolong(mhuc, mhuf) - - mhrf = assemble(conj(TestFunction(Vfine)) * dx) - mhrc = Cofunction(Vcoarse.dual()) - - tm.restrict(mhrf, mhrc) - - assert abs( - (assemble(action(mhrc, mhuc)) - assemble(action(mhrf, mhuf))) - / assemble(action(mhrf, mhuf)) - ) <= 1e-12 - assert abs( - (assemble(action(rc, u_coarse)) - assemble(action(mhrc, mhuc))) - / assemble(action(mhrc, mhuc)) - ) <= 1e-12 - - -@pytest.mark.parallel([1, 2]) -@pytest.mark.skipnetgen -def test_restrict_CG1(amh, atm): # pylint: disable=W0621 - """ - Test restriction with CG1 - """ - V_coarse = FunctionSpace(amh[0], "CG", 1) - V_fine = FunctionSpace(amh[-1], "CG", 1) +@pytest.mark.parallel([1, 2, 4]) +def test_restrict_CG1(mh): + """Test restriction with CG1""" + V_coarse = FunctionSpace(mh[0], "CG", 1) + V_fine = FunctionSpace(mh[-1], "CG", 1) u_coarse = Function(V_coarse) u_fine = Function(V_fine) xc, *_ = SpatialCoordinate(V_coarse.mesh()) u_coarse.interpolate(xc) - atm.prolong(u_coarse, u_fine) + prolong(u_coarse, u_fine) rf = assemble(conj(TestFunction(V_fine)) * dx) rc = Cofunction(V_coarse.dual()) - atm.restrict(rf, rc) + restrict(rf, rc) assert np.allclose( assemble(action(rc, u_coarse)), @@ -222,24 +372,21 @@ def test_restrict_CG1(amh, atm): # pylint: disable=W0621 ) -@pytest.mark.parallel([1, 2]) -@pytest.mark.skipnetgen -def test_restrict_DG0(amh, atm): # pylint: disable=W0621 - """ - Test restriction with DG0 - """ - V_coarse = FunctionSpace(amh[0], "DG", 0) - V_fine = FunctionSpace(amh[-1], "DG", 0) +@pytest.mark.parallel([1, 2, 4]) +def test_restrict_DG0(mh): + """Test restriction with DG0""" + V_coarse = FunctionSpace(mh[0], "DG", 0) + V_fine = FunctionSpace(mh[-1], "DG", 0) u_coarse = Function(V_coarse) u_fine = Function(V_fine) xc, *_ = SpatialCoordinate(V_coarse.mesh()) u_coarse.interpolate(xc) - atm.prolong(u_coarse, u_fine) + prolong(u_coarse, u_fine) rf = assemble(conj(TestFunction(V_fine)) * dx) rc = Cofunction(V_coarse.dual()) - atm.restrict(rf, rc) + restrict(rf, rc) assert np.allclose( assemble(action(rc, u_coarse)), @@ -249,13 +396,10 @@ def test_restrict_DG0(amh, atm): # pylint: disable=W0621 @pytest.mark.parallel([1, 2]) -@pytest.mark.skipnetgen -def test_mg_jacobi(amh, atm): # pylint: disable=W0621 - """ - Test multigrid with jacobi smoothers - """ - V = FunctionSpace(amh[-1], "CG", 1) - x = SpatialCoordinate(amh[-1]) +def test_mg_jacobi(mh): + """Test multigrid with jacobi smoothers""" + V = FunctionSpace(mh[-1], "CG", 1) + x = SpatialCoordinate(mh[-1]) u_ex = Function(V).interpolate(sin(2 * pi * x[0]) * sin(2 * pi * x[1])) u = Function(V) v = TestFunction(V) @@ -280,19 +424,15 @@ def test_mg_jacobi(amh, atm): # pylint: disable=W0621 problem = NonlinearVariationalProblem(F, u, bc) solver = NonlinearVariationalSolver(problem, solver_parameters=params) - solver.set_transfer_manager(atm) solver.solve() assert errornorm(u_ex, u) <= 1e-8 @pytest.mark.parallel([1, 2]) -@pytest.mark.skipnetgen -@pytest.mark.parametrize("params", ["jacobi", "asm", "patch"]) -def test_mg_patch(amh, atm, params): # pylint: disable=W0621 - """ - Test multigrid with patch relaxation - """ - if params == "jacobi": +@pytest.mark.parametrize("backend", ["jacobi", "patch", "tinyasm"]) +def test_mg_patch(mh, backend): + """Test multigrid with patch relaxation""" + if backend == "jacobi": solver_params = { "mat_type": "matfree", "ksp_type": "cg", @@ -308,7 +448,7 @@ def test_mg_patch(amh, atm, params): # pylint: disable=W0621 "pc_type": "lu", }, } - elif params == "patch": + elif backend == "patch": solver_params = { "mat_type": "matfree", "ksp_type": "cg", @@ -347,13 +487,14 @@ def test_mg_patch(amh, atm, params): # pylint: disable=W0621 "ksp_max_it": 1, "pc_type": "python", "pc_python_type": "firedrake.ASMStarPC", - "pc_star_backend": "tinyasm", + "pc_star_backend": backend, }, "mg_coarse": {"ksp_type": "preonly", "pc_type": "lu"}, } - V = FunctionSpace(amh[-1], "CG", 1) - x = SpatialCoordinate(amh[-1]) + mesh = mh[-1] + V = FunctionSpace(mesh, "CG", 1) + x = SpatialCoordinate(mesh) u_ex = Function(V).interpolate(sin(2 * pi * x[0]) * sin(2 * pi * x[1])) u = Function(V) v = TestFunction(V) @@ -363,7 +504,20 @@ def test_mg_patch(amh, atm, params): # pylint: disable=W0621 problem = NonlinearVariationalProblem(F, u, bc) solver = NonlinearVariationalSolver(problem, solver_parameters=solver_params) - solver.set_transfer_manager(atm) - solver.solve() + pc = solver.snes.ksp.pc + assert pc.getType() == "mg" + assert pc.getMGLevels() == len(mh) assert errornorm(u_ex, u) <= 1e-8 + + +def test_deprecated_adaptive_aliases(): + """The deprecated aliases warn, and forward their arguments.""" + mesh = UnitSquareMesh(2, 2) + with pytest.warns(FutureWarning): + mh = AdaptiveMeshHierarchy(mesh, nested=False) + assert mh[-1] is mesh + assert not mh.nested + + with pytest.warns(FutureWarning): + assert isinstance(AdaptiveTransferManager(), TransferManager) diff --git a/tests/firedrake/multigrid/test_snes_adapt.py b/tests/firedrake/multigrid/test_snes_adapt.py index 8fef8742d1..7819463e50 100644 --- a/tests/firedrake/multigrid/test_snes_adapt.py +++ b/tests/firedrake/multigrid/test_snes_adapt.py @@ -42,7 +42,6 @@ def mark_cells(ctx, current_solution): v = TestFunction(V) problem = NonlinearVariationalProblem((u - 1.0)*v*dx, u) solver = NonlinearVariationalSolver(problem, marking_callback=mark_cells) - solver.set_transfer_manager(AdaptiveTransferManager()) dm = solver.snes.getDM() with dmhooks.add_hooks(dm, solver, appctx=solver._ctx): @@ -62,14 +61,14 @@ def mark_cells(ctx, current_solution): @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 2]) def test_snes_adapt_sequence_with_adaptive_multigrid(): from netgen.occ import WorkPlane, Axes, OCCGeometry, X, Z rect1 = WorkPlane(Axes((0, 0, 0), n=Z, h=X)).Rectangle(1, 2).Face() rect2 = WorkPlane(Axes((0, 1, 0), n=Z, h=X)).Rectangle(2, 1).Face() mesh = Mesh(OCCGeometry(rect1 + rect2, dim=2).GenerateMesh(maxh=0.8)) - amh = AdaptiveMeshHierarchy(mesh) - atm = AdaptiveTransferManager() + mh = MeshHierarchy(mesh) V = FunctionSpace(mesh, "CG", 1) old_dim = V.dim() @@ -133,16 +132,15 @@ def mark_cells(ctx, current_solution): solver = LinearVariationalSolver(problem, solver_parameters=params, marking_callback=mark_cells) - solver.set_transfer_manager(atm) u_adapted = solver.solve() adapted_mesh = u_adapted.function_space().mesh() hierarchy, level = get_level(adapted_mesh) assert seen[0] == mesh - assert hierarchy is amh + assert hierarchy is mh assert level == refinements - assert len(amh) == refinements + 1 + assert len(mh) == refinements + 1 assert adapted_mesh is not mesh assert u_adapted is not uh assert u_adapted.function_space().dim() > old_dim