diff --git a/docs/beginner/create_xdmf.md b/docs/beginner/create_xdmf.md index b01f10993..18112ec74 100644 --- a/docs/beginner/create_xdmf.md +++ b/docs/beginner/create_xdmf.md @@ -122,8 +122,9 @@ Then: ### How vertex data is obtained -The compatibility group writer uses PETSc's `createInterpolation` to project -data to mesh vertices: +The compatibility group writer projects data to mesh vertices with the +interpolation matrix PETSc's `createInterpolation` builds (UW3 assembles its own +copy of it, which also works on a rank that owns none of its cells' nodes): - **P1 continuous variables**: DOFs match mesh vertices exactly, so the owned global-vector data is written directly to `/vertex_fields/` — no interpolation needed. diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index fe3776599..39f6ee644 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -181,7 +181,12 @@ path uses**: 2. a **mesh-owned** coarse tail (`mesh._custom_mg_coarse_meshes`, what `mesh.adapt()` leaves on a refinement child) is picked up opportunistically, with the barycentric→RBF builder fallback, so a failed build degrades to GAMG - rather than crashing the solve. + rather than crashing the solve. In parallel the retry and the degrade happen + on a verdict every rank holds: a builder failure, which is agreed across + ranks as soon as the builder returns, or a check on allreduced sizes + (`custom_mg.CollectiveBuildError`). Any other exception is one rank's while + the others are inside the build's collectives, so it is raised rather than + retried alone (#886). Step 2 used to be unreachable from here. The standard path's injection hook runs *after* the rotated dispatch has already returned, and the rotated builder diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index 74075e043..fbaf93892 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/src/underworld3/cython/petsc_compat.h @@ -366,3 +366,245 @@ PetscErrorCode UW_DMCreateBdIntegralSandbox(DM src, DM *sandbox) *sandbox = sdm; PetscFunctionReturn(PETSC_SUCCESS); } + +// Interpolation between two DMs on ONE plex (a coordinate DM and a basis DM +// cloned from it, a variable's sub-DM and a scratch clone of the mesh): the +// matrix DMCreateInterpolation builds through DMPlexComputeInterpolatorNested +// (isRefined = PETSC_FALSE), copied from PETSc v3.25.0 with one change (#902), +// under PETSc's licence: +// +// Copyright (c) 1991-2025, UChicago Argonne, LLC and the PETSc Developers and +// Contributors +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// The change: PETSc runs its preallocation loop only on a rank with local +// rows or columns (`if (locRows || locCols)`, since b4ec2976050, v3.24.0) but +// its fill loop on every rank that holds cells. A shared point belongs to the +// highest rank holding it, so a rank whose cells lie inside higher ranks' +// regions owns no rows and no columns; its cells' entries then reach owners +// that never preallocated them ("New nonzero caused a malloc" in +// MatAssemblyEnd) and the ranks that did not raise wait in the assembly's +// collective. Here both loops run over the same cells. Finite-element fields +// only, which are the only ones UW3 interpolates this way, and one plex only: +// PETSc dispatches here on the two DMs sharing their plex, and a pair on two +// plexes is refused on every rank. Drop this for DMCreateInterpolation once +// PETSc runs the two loops alike. +PetscErrorCode UW_DMCreateInterpolationWithinPlex(DM dmc, DM dmf, Mat *In) +{ + PetscSection fsection, fglobalSection, csection, cglobalSection; + PetscScalar *elemMat, *vals; + PetscFE *fe; + PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; + PetscInt cTotDim = 0, rTotDim = 0, m, n, locRows, locCols; + MatType mtype; + PetscBool ismatis; + Mat preallocator; + + PetscFunctionBeginUser; + { + // clones share their plex, cone array included; the verdict is + // agreed, since a rank holding no points cannot tell two plexes apart + PetscInt *conesC, *conesF, other; + + PetscCall(DMPlexGetCones(dmc, &conesC)); + PetscCall(DMPlexGetCones(dmf, &conesF)); + other = conesC != conesF; + PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &other, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dmc))); + PetscCheck(!other, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_INCOMP, + "UW_DMCreateInterpolationWithinPlex: the two DMs are not on one plex"); + } + PetscCall(DMGetGlobalSection(dmf, &fglobalSection)); + PetscCall(PetscSectionGetConstrainedStorageSize(fglobalSection, &m)); + PetscCall(DMGetGlobalSection(dmc, &cglobalSection)); + PetscCall(PetscSectionGetConstrainedStorageSize(cglobalSection, &n)); + PetscCall(DMGetMatType(dmc, &mtype)); + PetscCall(PetscStrcmp(mtype, MATIS, &ismatis)); + PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), In)); + PetscCall(MatSetSizes(*In, m, n, PETSC_DETERMINE, PETSC_DETERMINE)); + PetscCall(MatSetType(*In, ismatis ? MATAIJ : mtype)); + + PetscCall(DMGetDimension(dmf, &dim)); + PetscCall(DMGetLocalSection(dmf, &fsection)); + PetscCall(DMGetLocalSection(dmc, &csection)); + PetscCall(PetscSectionGetNumFields(fsection, &Nf)); + PetscCall(DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd)); + PetscCall(PetscCalloc1(Nf, &fe)); + for (f = 0; f < Nf; ++f) { + PetscObject obj, objc; + PetscClassId id, idc; + PetscInt rNb, cNb; + + PetscCall(DMGetField(dmf, f, NULL, &obj)); + PetscCall(PetscObjectGetClassId(obj, &id)); + PetscCall(DMGetField(dmc, f, NULL, &objc)); + PetscCall(PetscObjectGetClassId(objc, &idc)); + PetscCheck(id == PETSCFE_CLASSID && idc == PETSCFE_CLASSID, PetscObjectComm((PetscObject)dmc), + PETSC_ERR_SUP, "UW_DMCreateInterpolationWithinPlex: field %" PetscInt_FMT " is not a PetscFE", f); + fe[f] = (PetscFE)obj; + PetscCall(PetscFEGetDimension(fe[f], &rNb)); + PetscCall(PetscFEGetDimension((PetscFE)objc, &cNb)); + rTotDim += rNb; + cTotDim += cNb; + } + PetscCall(PetscCalloc1(rTotDim * cTotDim, &elemMat)); + for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { + PetscDualSpace Qref; + PetscQuadrature quad; + const PetscReal *qpoints, *qweights; + PetscReal *points; + PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; + + // the fine field's dual basis: every functional's quadrature points + PetscCall(PetscFEGetDualSpace(fe[fieldI], &Qref)); + PetscCall(PetscFEGetNumComponents(fe[fieldI], &Nc)); + PetscCall(PetscDualSpaceGetDimension(Qref, &fpdim)); + for (i = 0; i < fpdim; ++i) { + PetscCall(PetscDualSpaceGetFunctional(Qref, i, &quad)); + PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Np, NULL, NULL)); + npoints += Np; + } + PetscCall(PetscMalloc1(npoints * dim, &points)); + for (i = 0, k = 0; i < fpdim; ++i) { + PetscCall(PetscDualSpaceGetFunctional(Qref, i, &quad)); + PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Np, &qpoints, NULL)); + for (p = 0; p < Np; ++p, ++k) + for (d = 0; d < dim; ++d) points[k * dim + d] = qpoints[p * dim + d]; + } + for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { + PetscObject obj; + PetscFE fec; + PetscTabulation T = NULL; + PetscInt NcJ, cpdim, j, qNc; + + PetscCall(DMGetField(dmc, fieldJ, NULL, &obj)); + fec = (PetscFE)obj; + PetscCall(PetscFEGetNumComponents(fec, &NcJ)); + PetscCall(PetscFEGetDimension(fec, &cpdim)); + // fields only interpolate themselves: the coarse basis evaluated at + // the fine functionals' points, weighted by the functionals + if (fieldI == fieldJ) { + PetscCheck(Nc == NcJ, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %" PetscInt_FMT " does not match coarse field %" PetscInt_FMT, Nc, NcJ); + PetscCall(PetscFECreateTabulation(fec, 1, npoints, points, 0, &T)); + for (i = 0, k = 0; i < fpdim; ++i) { + PetscCall(PetscDualSpaceGetFunctional(Qref, i, &quad)); + PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Np, NULL, &qweights)); + PetscCheck(qNc == NcJ, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %" PetscInt_FMT " does not match coarse field %" PetscInt_FMT, qNc, NcJ); + for (p = 0; p < Np; ++p, ++k) + for (j = 0; j < cpdim; ++j) + for (c = 0; c < Nc; ++c) elemMat[(offsetI + i) * cTotDim + offsetJ + j] += T->T[0][k * cpdim * NcJ + j * Nc + c] * qweights[p * qNc + c]; + } + PetscCall(PetscTabulationDestroy(&T)); + } + offsetJ += cpdim; + } + offsetI += fpdim; + PetscCall(PetscFree(points)); + } + + // Preallocate from every cell this rank holds, as the fill below does, + // whether or not the rank owns any rows or columns (the one change). + PetscCall(MatGetLocalSize(*In, &locRows, &locCols)); + PetscCall(MatCreate(PetscObjectComm((PetscObject)*In), &preallocator)); + PetscCall(MatSetType(preallocator, MATPREALLOCATOR)); + PetscCall(MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE)); + PetscCall(MatSetUp(preallocator)); + PetscCall(PetscCalloc1(rTotDim * cTotDim, &vals)); + for (c = cStart; c < cEnd; ++c) + PetscCall(DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, PETSC_FALSE, dmc, csection, cglobalSection, PETSC_FALSE, preallocator, c, vals, INSERT_VALUES)); + PetscCall(PetscFree(vals)); + PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY)); + PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY)); + PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *In)); + PetscCall(MatDestroy(&preallocator)); + + PetscCall(MatZeroEntries(*In)); + for (c = cStart; c < cEnd; ++c) + PetscCall(DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, PETSC_FALSE, dmc, csection, cglobalSection, PETSC_FALSE, *In, c, elemMat, INSERT_VALUES)); + PetscCall(PetscFree(fe)); + PetscCall(PetscFree(elemMat)); + PetscCall(MatAssemblyBegin(*In, MAT_FINAL_ASSEMBLY)); + PetscCall(MatAssemblyEnd(*In, MAT_FINAL_ASSEMBLY)); + PetscCall(MatViewFromOptions(*In, NULL, "-interp_mat_view")); + PetscFunctionReturn(PETSC_SUCCESS); +} + +// The coordinate space of a plex at `degree`, as +// DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_TRUE) builds it, +// with the coordinates carried across by UW_DMCreateInterpolationWithinPlex. +// PETSc's projection (DMSetCoordinateDisc) interpolates through +// DMCreateInterpolation between the old and the new coordinate DM, two DMs on +// this plex, so a degree >= 2 space could not be built on a rank that holds +// cells but owns none of their nodes (#902). As in PETSc: the same space is a +// copy; coordinates injected into a higher coordinate dimension are not +// projected; a cell type PETSc builds no Lagrange space for (tensor prisms) +// keeps its coordinates. Isoperiodic coordinates keep PETSc's own projection, +// which is local there, without the matrix; maxCell/L periodicity projects +// the non-localized coordinates through the matrix, as here. +PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(DM dm, PetscInt degree) +{ + PetscInt nfaceSF = 0; + DM cdmOld, cdmNew; + Vec xOld, xNew, x; + PetscFE discOld, disc; + PetscDualSpace dsOld, ds; + PetscBool same; + + PetscFunctionBeginUser; + PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &nfaceSF, NULL)); + if (nfaceSF) { + PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_TRUE)); + PetscFunctionReturn(PETSC_SUCCESS); + } + PetscCall(DMGetCoordinateDM(dm, &cdmOld)); + PetscCall(PetscObjectReference((PetscObject)cdmOld)); + PetscCall(DMGetCoordinates(dm, &xOld)); + PetscCall(PetscObjectReference((PetscObject)xOld)); + // installs the new coordinate DM and leaves the coordinates to us + PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_FALSE, PETSC_FALSE)); + PetscCall(DMGetCoordinates(dm, &x)); + PetscCall(DMGetCoordinateDM(dm, &cdmNew)); + // no new coordinate DM: PETSc installed no space for this cell type + if (x == xOld && cdmNew != cdmOld) { + PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld)); + PetscCall(DMGetField(cdmNew, 0, NULL, (PetscObject *)&disc)); + PetscCall(PetscFEGetDualSpace(discOld, &dsOld)); + PetscCall(PetscFEGetDualSpace(disc, &ds)); + PetscCall(PetscDualSpaceEqual(dsOld, ds, &same)); + PetscCall(DMCreateGlobalVector(cdmNew, &xNew)); + if (same) { + PetscCall(VecCopy(xOld, xNew)); + } else { + Mat In; + + PetscCall(UW_DMCreateInterpolationWithinPlex(cdmOld, cdmNew, &In)); + PetscCall(MatMult(In, xOld, xNew)); + PetscCall(MatDestroy(&In)); + } + PetscCall(DMSetCoordinates(dm, xNew)); + PetscCall(VecDestroy(&xNew)); + } + PetscCall(VecDestroy(&xOld)); + PetscCall(DMDestroy(&cdmOld)); + PetscFunctionReturn(PETSC_SUCCESS); +} diff --git a/src/underworld3/cython/petsc_discretisation.pyx b/src/underworld3/cython/petsc_discretisation.pyx index a15613057..364be1ecb 100644 --- a/src/underworld3/cython/petsc_discretisation.pyx +++ b/src/underworld3/cython/petsc_discretisation.pyx @@ -181,6 +181,41 @@ def petsc_dm_set_regular_refinement(dm, regular=True): return +def petsc_dm_create_interpolation_within_plex(dmc, dmf): + """The interpolation matrix from the field of ``dmc`` to the field of + ``dmf``, two DMs on ONE plex (clones or sub-DMs of one mesh DM). + + It is the matrix ``dmc.createInterpolation(dmf)`` returns, assembled + also on a rank that holds cells but owns none of their degrees of + freedom, where PETSc's own preallocation skips the rank and the + assembly fails on some ranks and hangs the others (#902; see + ``UW_DMCreateInterpolationWithinPlex`` in ``petsc_compat.h``). + Collective. + """ + cdef DM c_dmc = dmc + cdef DM c_dmf = dmf + cdef Mat interp = PETSc.Mat() + CHKERRQ( UW_DMCreateInterpolationWithinPlex(c_dmc.dm, c_dmf.dm, &interp.mat) ) + return interp + + +def petsc_dm_create_coordinate_space_within_plex(dm, degree): + """Give ``dm`` a Lagrange coordinate space of ``degree``, its coordinates + carried across from the current space. + + It is ``dm.createCoordinateSpace(degree, False, True)``, with the + coordinates interpolated by UW3's copy of PETSc's interpolator, so a + degree >= 2 space is also built on a rank that holds cells but owns + none of their nodes (#902; see + ``UW_DMPlexCreateCoordinateSpaceWithinPlex`` in ``petsc_compat.h``). + Isoperiodic coordinates are the exception: PETSc projects them + locally, with no matrix, and they keep its route. Collective. + """ + cdef DM c_dm = dm + cdef PetscInt c_degree = degree + CHKERRQ( UW_DMPlexCreateCoordinateSpaceWithinPlex(c_dm.dm, c_degree) ) + + # This is not cython, does it need to be here or in discretisation.py ? diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index 4c827d2e2..337f157e3 100644 --- a/src/underworld3/cython/petsc_extras.pxi +++ b/src/underworld3/cython/petsc_extras.pxi @@ -58,6 +58,8 @@ cdef extern from "petsc_compat.h": PetscErrorCode UW_DMForceCoordinateField(PetscDM) PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*) PetscErrorCode UW_DMCreateBdIntegralSandbox(PetscDM, PetscDM*) + PetscErrorCode UW_DMCreateInterpolationWithinPlex(PetscDM, PetscDM, PetscMat*) + PetscErrorCode UW_DMPlexCreateCoordinateSpaceWithinPlex(PetscDM, PetscInt) cdef extern from "petsc.h" nogil: PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[]) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 453781317..f903a555a 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2754,8 +2754,14 @@ def nuke_coords_and_rebuild( # dual space without proper point subspaces, causing # DMPlexComputeBdIntegral to segfault/deadlock (issue #96). # createCoordinateSpace builds the FE internally with correct - # subspace initialisation. - self.dm.createCoordinateSpace(self.degree, False, True) + # subspace initialisation. The coordinates are carried across by + # UW3's copy of PETSc's interpolator: PETSc's own projection of a + # degree >= 2 space cannot be assembled on a rank that holds cells + # but owns none of their nodes (#902). + from underworld3.cython.petsc_discretisation import ( + petsc_dm_create_coordinate_space_within_plex, + ) + petsc_dm_create_coordinate_space_within_plex(self.dm, self.degree) # Issue #96 fix: Force coordinate field creation and strip # boundary labels from the coordinate DM. createCoordinateSpace @@ -5729,24 +5735,29 @@ def _cell_node_indices(self, degree, continuous): return self._cell_node_array[key] def _get_coords_for_basis(self, degree, continuous): + """The coordinates of the nodes of a degree-``degree`` Lagrange field, an + ``(n, cdim)`` array in the local layout of :meth:`_basis_coordinate_dm` + (ghost nodes included). + + The interpolation matrix comes from UW3's copy of PETSc's nested + interpolator, which assembles also on a rank that holds cells but owns + none of their nodes (#902; ``petsc_dm_create_interpolation_within_plex``). """ - This function returns the vertex array for the - provided variable. If the array does not already exist, - it is first created and then returned. - """ + from underworld3.cython.petsc_discretisation import ( + petsc_dm_create_interpolation_within_plex, + ) dmold = self.dm.getCoordinateDM() dmnew = self._basis_coordinate_dm(degree, continuous) - matInterp, vecScale = dmold.createInterpolation(dmnew) + matInterp = petsc_dm_create_interpolation_within_plex(dmold, dmnew) coordsOld = self.dm.getCoordinates() coordsNewL = dmnew.getLocalVec() coordsNewG = dmnew.getGlobalVec() matInterp.mult(coordsOld, coordsNewG) dmnew.globalToLocal(coordsNewG, coordsNewL) - arr = coordsNewL.array - arrcopy = arr.reshape(-1, self.cdim).copy() + arrcopy = coordsNewL.array.reshape(-1, self.cdim).copy() dmnew.restoreGlobalVec(coordsNewG) dmnew.restoreLocalVec(coordsNewL) @@ -5754,8 +5765,6 @@ def _get_coords_for_basis(self, degree, continuous): # they accumulate until Python GC runs — noticeable in long adapt # loops that re-fill the coord cache per variable. matInterp.destroy() - if vecScale is not None: - vecScale.destroy() dmnew.destroy() return arrcopy diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index d8922776e..b0551e484 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -1,8 +1,9 @@ """Project or prolong MeshVariable data to a different polynomial degree. -Uses a scratch PETSc DM with ``createInterpolation`` — no MeshVariable -is created, the mesh DM is never modified, and all scratch objects are -destroyed before returning. +Uses a scratch PETSc DM and the interpolation matrix between it and the +variable's sub-DM (UW3's copy of PETSc's ``createInterpolation``, #902) — no +MeshVariable is created, the mesh DM is never modified, and all scratch +objects are destroyed before returning. Typical use cases ----------------- @@ -81,9 +82,19 @@ def project_to_degree( iset, subdm_src = mesh.dm.createSubDM(mesh_var.field_id) # --- interpolation matrix --- - # PETSc: source_subdm.createInterpolation(target_dm) returns a matrix - # whose mult() maps source global vec → target global vec. - interp, _scale = subdm_src.createInterpolation(dm_scratch) + # Its mult() maps source global vec → target global vec. UW3's copy of + # PETSc's nested interpolator: PETSc's own (createInterpolation) cannot be + # assembled on a rank that holds cells but owns none of their DOFs (#902). + from underworld3.cython.petsc_discretisation import ( + petsc_dm_create_interpolation_within_plex, + ) + interp = petsc_dm_create_interpolation_within_plex(subdm_src, dm_scratch) + + # TODO(BUG): a discontinuous source projected to a continuous target sums + # the adjacent cells' values at each shared node (2-7x at P1 vertices): the + # matrix inserts one entry per cell into different source columns. The + # XDMF writer is not affected (it projects a discontinuous field to DG0). + # See issue #905. # --- apply --- g_src = subdm_src.createGlobalVec() @@ -103,8 +114,6 @@ def project_to_degree( # --- cleanup --- for obj in (interp, g_src, g_dst, dm_scratch, fe_target, subdm_src, iset): obj.destroy() - if _scale is not None: - _scale.destroy() return result diff --git a/src/underworld3/mpi.py b/src/underworld3/mpi.py index 97761d3d2..e804375b0 100644 --- a/src/underworld3/mpi.py +++ b/src/underworld3/mpi.py @@ -55,6 +55,14 @@ size = comm.size rank = comm.rank +# The communicator of collective error verdicts, and of nothing else. A verdict +# exchanged on ``comm`` can pair with a collective of the same kind that +# another rank is still inside (a rank that raised before it, the others in +# it): the two then trade unrelated data, and the ranks leave with different +# results. On its own communicator a verdict meets only verdicts, and a rank +# that never reaches it leaves the others waiting at a known place (#886). +verdict_comm = comm.Dup() + # State tracking for selective execution _in_selective_ranks = False _selective_executing_ranks = None diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 6e5d707b2..ec8654ac1 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -154,7 +154,12 @@ def _trace_interior_coords(solver, degree): mesh.dim, cdim, mesh.isSimplex, mesh.qdegree, prefix, PETSc.COMM_SELF) dmnew.setField(0, fe) dmnew.createDS() - mat_interp, vec_scale = dmold.createInterpolation(dmnew) + # UW3's copy of PETSc's nested interpolator: PETSc's own cannot be + # assembled on a rank that holds cells but owns none of their nodes (#902) + from underworld3.cython.petsc_discretisation import ( + petsc_dm_create_interpolation_within_plex, + ) + mat_interp = petsc_dm_create_interpolation_within_plex(dmold, dmnew) coords_new_g = dmnew.getGlobalVec() coords_new_l = dmnew.getLocalVec() mat_interp.mult(solver.dm.getCoordinates(), coords_new_g) @@ -171,8 +176,6 @@ def _trace_interior_coords(solver, degree): dmnew.restoreGlobalVec(coords_new_g) dmnew.restoreLocalVec(coords_new_l) mat_interp.destroy() - if vec_scale is not None: - vec_scale.destroy() fe.destroy() dmnew.destroy() return out diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index dd5081cd6..3b49d6362 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -224,6 +224,19 @@ def _drop_structural_zeros(P_csr, tol=1e-12): _BUILDERS = {"barycentric": barycentric_prolongation, "rbf": rbf_prolongation} +class CollectiveBuildError(RuntimeError): + """A hierarchy-build failure raised on every rank together, from a verdict + every rank holds (a builder that failed on any rank, a check on allreduced + sizes), so the ranks can retry or degrade together (#886). + ``builder_failure`` says the transfer builder itself failed, the one case + a different builder can change; any other check failing again under + another builder, the build degrades at once.""" + + def __init__(self, message, builder_failure=False): + super().__init__(message) + self.builder_failure = builder_failure + + def _require_serial(where): """Custom-P transfers are serial-only (experimental). The reduced maps use rank-local DOF indices and the prolongations assemble as serial AIJ; at np>1 @@ -1107,6 +1120,25 @@ def _node_transfer_no_repair(builder, coarse_coords, fine_coords): return builder(coarse_coords, fine_coords) +def _agreed_node_transfer(builder, coarse_coords, fine_coords, comm): + """The rank-local node transfer of a parallel build, with one verdict for + every rank: if the builder raised on any rank, :class:`CollectiveBuildError` + is raised on all of them. The builders do no communication, so every rank + reaches the verdict; a rank that kept its own error would go on to retry + while the others assemble the transfer, and the two would meet in different + collectives (#886).""" + problem = None + try: + Pn = _drop_structural_zeros( + _node_transfer_no_repair(builder, coarse_coords, fine_coords)) + except Exception as exc: # sanctioned: the verdict below raises it on every rank + Pn, problem = None, f"rank {comm.getRank()}: {exc}" + problems = [p for p in comm.tompi4py().allgather(problem) if p is not None] + if problems: + raise CollectiveBuildError("; ".join(problems), builder_failure=True) + return Pn + + def _assemble_node_transfer(Pn, coarse_cols, fine_layout, coarse_layout, ncomp, comm): """A node-level transfer ``Pn`` (scipy CSR, fine local nodes by coarse @@ -1156,8 +1188,7 @@ def _build_parallel_transfer(coarse_coords, fine_coords, coarse_layout, coarse ``l2g`` to global columns (off-rank columns are fine for MPIAIJ). Constrained coarse DOFs (``l2g == -1``) drop out -> reduced->reduced.""" l2g_c = coarse_layout.l2g - Pn = _drop_structural_zeros( - _node_transfer_no_repair(builder, coarse_coords, fine_coords)) # (n_f_nodes, n_c_nodes), local + Pn = _agreed_node_transfer(builder, coarse_coords, fine_coords, comm) # (n_f_nodes, n_c_nodes), local return _assemble_node_transfer( Pn, np.asarray(l2g_c, dtype=np.int64).reshape(-1, ncomp), fine_layout, coarse_layout, ncomp, comm) @@ -1199,8 +1230,7 @@ def _build_crosspart_transfer(coarse_coords, fine_coords, coarse_layout, coords_u, cols_u = _gather_coarse_cloud(coarse_coords, coarse_layout, ncomp, comm) - Pn = _drop_structural_zeros( - _node_transfer_no_repair(builder, coords_u, fine_coords)) # (n_f_local, Nu) + Pn = _agreed_node_transfer(builder, coords_u, fine_coords, comm) # (n_f_local, Nu) return _assemble_node_transfer(Pn, np.asarray(cols_u, dtype=np.int64), fine_layout, coarse_layout, ncomp, comm) @@ -1743,14 +1773,15 @@ def build(self, solver): level_comm = solver.dm.comm.tompi4py() nodes, dofs = level_comm.allreduce(nodes), level_comm.allreduce(dofs) nc = dofs // nodes - if dofs % nodes != 0: - raise RuntimeError( + if dofs % nodes != 0: # allreduced: every rank raises + raise CollectiveBuildError( f"level {k}: full DOFs {dofs} not divisible by nodes {nodes}") coords.append(c) ncomp.append(nc) - if len(set(ncomp)) != 1: - raise RuntimeError(f"inconsistent component counts across levels: {ncomp}") + if len(set(ncomp)) != 1: # from allreduced counts: every rank raises + raise CollectiveBuildError( + f"inconsistent component counts across levels: {ncomp}") nc = ncomp[0] # Operator-faithful check: the finest reduced map must span exactly the @@ -2087,8 +2118,8 @@ def _assert_finest_matches_operator(solver, finest_map, parallel): int(finest_map.rend - finest_map.rstart))) else: red_n = int(len(finest_map)) # r2f length = reduced global size - if red_n != op_n: - raise RuntimeError( + if red_n != op_n: # allreduced: every rank raises + raise CollectiveBuildError( f"custom_mg: finest reduced-map size {red_n} != assembled operator " f"size {op_n}. The DM global section disagrees with the operator " f"(an adapt-child section inconsistency); the finest transfer would " @@ -2163,12 +2194,12 @@ def __init__(self, dm, fine_mesh): self.qdegree = fine_mesh.qdegree # The gmsh-imported BASE level carries section-only coordinates (its - # coordinate field is a PetscContainer, no PetscFE), and - # DMCreateInterpolation from such a source silently returns a ZERO - # matrix (measured) — _get_coords_for_basis would then hand every - # level-0 node the coordinate (0,0) and the transfer build collapses. - # refine() gives the child levels an FE coordinate space; give the - # clone's base the same: P1 Lagrange on the identical vertex layout. + # coordinate field is a PetscContainer, no PetscFE). PETSc's + # DMCreateInterpolation from such a source silently returned a ZERO + # matrix (measured), and UW3's copy, which _get_coords_for_basis now + # uses (#902), refuses a field that is not a PetscFE. refine() gives + # the child levels an FE coordinate space; give the clone's base the + # same: P1 Lagrange on the identical vertex layout. cdm = self.dm.getCoordinateDM() field = cdm.getField(0) fobj = field[0] if isinstance(field, tuple) else field @@ -2259,7 +2290,9 @@ def build_transfers(solver, field_id=None): is a DEMAND (build errors raise — the user registered it explicitly); the other two are PREFERENCES, built through the same opportunistic arm (barycentric, RBF retry, degrade to the solver's default preconditioner — - every step recorded in ``solver.pc_fallbacks``). + every step recorded in ``solver.pc_fallbacks``). In parallel only a builder + failure every rank raises together (:class:`CollectiveBuildError`) is + retried or degraded; any other failure is raised where it happens (#886). A refinement child carries ``mesh._custom_mg_coarse_meshes`` (the static coarse tail), so a :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]`` @@ -2370,6 +2403,13 @@ def build_transfers(solver, field_id=None): # barycentric fell back to GAMG at 23 its while RBF kept pc=mg at 2. # Falling back to GAMG loses the hierarchy entirely, so try the cheaper # degradation first. See #424. + # In parallel the retry and the degrade are collective decisions, so they + # are taken only on a verdict every rank holds (CollectiveBuildError: a + # builder failure, a check on allreduced sizes). Any other exception is + # one rank's: the others are inside the build's collectives, and a rank + # that caught it and retried would meet them in different ones (#886). + # It is raised. + _retryable = CollectiveBuildError if solver.dm.comm.getSize() > 1 else Exception _attempts = [builder] + (["rbf"] if builder != "rbf" else []) h = Ps = None for _i, _b in enumerate(_attempts): @@ -2382,9 +2422,14 @@ def build_transfers(solver, field_id=None): try: Ps = h.build(solver) break - except Exception as exc: # pragma: no cover - defensive + except _retryable as exc: # sanctioned: retried or degraded, and recorded import warnings - if _i + 1 < len(_attempts): + # another builder can only rescue a failure of the builder: an + # agreed check elsewhere (the provenance transfer, the operator + # sizes) would fail again after a full rebuild + retry = (exc.builder_failure if isinstance(exc, CollectiveBuildError) + else True) + if retry and _i + 1 < len(_attempts): solver._record_pc_fallback( "custom_mg.transfer_builder", requested=_b, @@ -2404,7 +2449,7 @@ def build_transfers(solver, field_id=None): continue solver._record_pc_fallback( "custom_mg.build", - requested=f"custom-P geometric MG ({' -> '.join(_attempts)})", + requested=f"custom-P geometric MG ({' -> '.join(_attempts[:_i + 1])})", installed="default preconditioner", reason="build_failed", detail=str(exc)) diff --git a/src/underworld3/utilities/fault_contact.py b/src/underworld3/utilities/fault_contact.py index f55501a68..71482b448 100644 --- a/src/underworld3/utilities/fault_contact.py +++ b/src/underworld3/utilities/fault_contact.py @@ -70,7 +70,9 @@ def _collective_raise(problem): where an error was meant. COLLECTIVE — every rank must reach it. """ if mpi.size > 1: - gathered = mpi.comm.allgather(problem) + # on the verdict communicator: a rank still inside a collective of + # the path cannot be paired with this exchange + gathered = mpi.verdict_comm.allgather(problem) problems = [p for p in gathered if p is not None] problem = problems[0] if problems else None if problem is not None: diff --git a/src/underworld3/utilities/fault_mg.py b/src/underworld3/utilities/fault_mg.py index a2f1b7008..65af717cc 100644 --- a/src/underworld3/utilities/fault_mg.py +++ b/src/underworld3/utilities/fault_mg.py @@ -461,9 +461,12 @@ def provenance_transfer_parallel(coarse, fine, pre, degree, continuous, nv = cv_c.shape[1] dim = Xc.shape[1] + # both checks below are on allreduced counts, so every rank raises + # together: build_transfers may retry or degrade on them (#886) + from underworld3.utilities.custom_mg import CollectiveBuildError owner = _fine_cell_of_node(fn, n_f) if m4.allreduce(int((owner < 0).sum())): - raise RuntimeError("fault_mg: a fine node belongs to no cell") + raise CollectiveBuildError("fault_mg: a fine node belongs to no cell") mine = np.flatnonzero(((l2g_f >= fstart) & (l2g_f < fend)).any(axis=1)) fcent = vxy_f[cv_f].mean(axis=1) if pre is fine: @@ -478,7 +481,7 @@ def provenance_transfer_parallel(coarse, fine, pre, degree, continuous, n_out = int(m4.allreduce(int((plam.min(axis=1) < -_INSIDE_TOL).sum()) if mine.size else 0)) if n_out: - raise RuntimeError( + raise CollectiveBuildError( f"fault_mg: {n_out} fine node(s) lie in no cell of the level as the " f"add-on partition left it; the surgery was not rank-local.") parent = np.asarray(pre._addon_parent, dtype=np.int32)[pcell] # (nm, 2) @@ -612,12 +615,16 @@ def split_tail_levels(levels, fault_specs, verbose=False): # partition ran once, for every level): no second partition split = add_fault(level, entries, verbose=verbose, partition=False, **kwargs) - except (ValueError, RuntimeError, NotImplementedError) as exc: - problem = str(exc) + except Exception as exc: # sanctioned: agreed below, the level stays welded + problem = f"{type(exc).__name__}: {exc}" # the split's refusals are collective; a failure on one rank after its # last collective is not, and every rank must stop at the same level or - # the ranks hold different hierarchies (a rank-local raise inside the - # split, before its own collectives, still hangs there: issue #874) + # the ranks hold different hierarchies. Any exception is caught, so an + # IndexError on one rank meets the same verdict (#886). A rank-local + # raise inside the split, before its own collectives, still hangs: the + # raising rank waits in the verdict, which has its own communicator so + # it cannot pair with the split's collectives, and the others wait in + # the split (issue #874) try: _collective_raise(problem) except RuntimeError as exc: diff --git a/tests/parallel/test_1089_fault_split_multigrid_parallel.py b/tests/parallel/test_1089_fault_split_multigrid_parallel.py index 2fd5489c7..187648001 100644 --- a/tests/parallel/test_1089_fault_split_multigrid_parallel.py +++ b/tests/parallel/test_1089_fault_split_multigrid_parallel.py @@ -45,6 +45,46 @@ def test_split_levels_converge_like_serial(): assert abs(slip - SERIAL_PEAK_SLIP) < 1.0e-4 * SERIAL_PEAK_SLIP, slip +def test_an_agreed_provenance_failure_degrades_every_rank(monkeypatch): + """The provenance transfer checks allreduced counts ("a fine node belongs + to no cell"), so every rank raises CollectiveBuildError together and the + multigrid build degrades on every rank, not on one (#886). A different + transfer builder cannot help, since the provenance transfer uses none, so + there is no retry with the RBF builder (which repeated the same failure + after a full rebuild). Here one fine node on rank 1 is made to belong to no + cell; the solve then runs on the default preconditioner, to the same slip.""" + from underworld3.utilities import fault_mg + original_owner = fault_mg._fine_cell_of_node + original_provenance = fault_mg.provenance_transfer_parallel + inside, calls = [False], [0] + + def counted(*args, **kwargs): + calls[0] += 1 + inside[0] = True + try: + return original_provenance(*args, **kwargs) + finally: + inside[0] = False + + def one_orphan_on_rank_one(fn, n_f): + owner = original_owner(fn, n_f) + if inside[0] and uw.mpi.rank == 1 and owner.size: + owner = owner.copy() + owner[0] = -1 + return owner + + monkeypatch.setattr(fault_mg, "_fine_cell_of_node", one_orphan_on_rank_one) + monkeypatch.setattr(fault_mg, "provenance_transfer_parallel", counted) + mesh = _faulted_box(2) + with pytest.warns(UserWarning, match="custom-P FMG build failed"): + stokes, info, slip = shear_fault_solve(mesh) + assert info["converged"] and info["velocity_pc"] == "GAMG", info + assert uw.mpi.comm.allgather(calls[0]) == [1] * uw.mpi.size + assert "custom_mg.build" in stokes.pc_fallbacks + assert "custom_mg.transfer_builder" not in stokes.pc_fallbacks + assert abs(slip - SERIAL_PEAK_SLIP) < 1.0e-4 * SERIAL_PEAK_SLIP, slip + + def test_the_welded_tail_runs_in_parallel_too(): info, slip = _solve(_faulted_box(2, split_levels=False)) assert info["velocity_pc"] == "custom-FMG" and info["converged"], info @@ -94,14 +134,17 @@ def centroid(q): assert np.abs(np.einsum("ij,ij->i", jumps, normals)).max() < 1.0e-8 -def test_a_failure_on_one_rank_stops_the_replay_on_every_rank(monkeypatch): +@pytest.mark.parametrize("error", [RuntimeError, IndexError]) +def test_a_failure_on_one_rank_stops_the_replay_on_every_rank(monkeypatch, error): """A failure of the split on one rank after its last collective (a final consistency check, rank-local) stopped the replay of the fault on that rank alone: the ranks kept different split levels and the next collective on the hierarchy met mismatched partners. The verdict is exchanged per level, so every rank keeps the same split levels and the levels from the failure down - stay welded. (A rank-local raise BEFORE the split's own collectives still - hangs inside the split; issue #874.)""" + stay welded. Any exception meets the verdict, not only the split's own + refusal types: an IndexError escaped it on one rank (#886). (A rank-local + raise BEFORE the split's own collectives still hangs, the raising rank in + the verdict and the others in the split; issue #874.)""" import underworld3.utilities.fault_split as fault_split original, calls = fault_split.add_fault, [0] @@ -109,7 +152,7 @@ def fails_on_one_rank_after_the_second_add_on_split(mesh, *args, **kwargs): calls[0] += 1 # 1: the finest; 2, 3: the add-on levels split = original(mesh, *args, **kwargs) if calls[0] == 3 and uw.mpi.rank == 1: - raise RuntimeError("an internal failure on one rank") + raise error("an internal failure on one rank") return split monkeypatch.setattr(fault_split, "add_fault", diff --git a/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py new file mode 100644 index 000000000..9ffe7d5a0 --- /dev/null +++ b/tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py @@ -0,0 +1,285 @@ +"""Interpolation within one mesh on a rank that holds cells but owns none of their DOFs. + +PETSc's ``DMPlexComputeInterpolatorNested`` builds the matrix between two DMs on +one plex (a coordinate DM and a basis DM, a variable's sub-DM and a scratch +clone). Since v3.24.0 it skips its preallocation pass on a rank with no local +rows and no local columns, but runs its fill pass over that rank's cells all +the same. A +shared point belongs to the highest rank that holds it, so a rank whose cells +sit wholly inside higher ranks' regions owns nothing. Its entries then reach +owners that never preallocated them: the owner raises "New nonzero caused a +malloc" inside ``MatAssemblyEnd``, and the rank that did not raise waits in the +assembly's own collective for ever. Found on kaiju at np 32 and 52 on the +124-cell coarsest level of the faulted annulus, where custom-P multigrid reads +each level's node coordinates (#902), and on a 3-D notch mesh at np 16. + +A shell partition gives rank 0 one cell of a 4 x 4 grid of triangles and the +other ranks the rest. An interior cell leaves rank 0 owning none of its +vertices and edges: PETSc's P2 matrices fail there, while its P1 matrices +escape, because every vertex pair of the cell is also an edge of another +rank's cell. A boundary cell leaves rank 0 owning its wall edge but none of +its vertices: there PETSc's P1 matrices fail, since the wall edge's vertex +pair occurs in no other cell, and its P2 matrices pass. At degree 3 rank 0 +owns its cell's interior node, so PETSc preallocates; those cases check the +copy's node placement. The tests check each case is set up, then +that the node coordinates (degree 1 to 3, the third with two nodes on each +edge) and the projection between degrees come back right on every rank, +through UW3's copy of the interpolator +(``petsc_dm_create_interpolation_within_plex``). A degree-2 mesh on the +partition pins the coordinate space PETSc projects in ``Mesh`` setup, and a +custom-P FMG solve whose coarse level is the grid, under its refinement, pins +the case inside a real multigrid build. + +Before the fix the ranks failed or hung inside PETSc. The thread timeout below +stops a rank waiting in Python; a rank held inside a PETSc collective keeps the +GIL, so a hang there is ended by the MPI supervisor's silence limit instead. + + mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1107_interpolation_rank_owning_nothing_parallel.py +""" +from enum import Enum + +import numpy as np +import pytest +from mpi4py import MPI +from petsc4py import PETSc + +import underworld3 as uw +from underworld3.function import project_to_degree +from underworld3.utilities import custom_mg + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, pytest.mark.tier_b, + pytest.mark.timeout(120, method="thread")] + +N = 4 +# Custom-P FMG iterations of the test below, at np 2, 3 and 4 (GAMG takes 12). A +# transfer built from coarse nodes displaced by 0.3 of a cell takes 5. +FMG_ITERATIONS = 4 +_WALLS = Enum("Boundaries", {"Walls": 1}) + + +def _cells_of(dm, v): + return {c for e in dm.getSupport(v) for c in dm.getSupport(e)} + + +def _owned_nothing_dm(where="interior"): + """The 4 x 4 grid of triangles, all on rank 0, partitioned so rank 0 keeps one + cell: an ``interior`` one (every edge shared with a cell of another rank), + or a ``boundary`` one with a single wall edge and every vertex also in + another cell. Its outer edges carry the label ``Walls``.""" + comm = PETSc.COMM_WORLD + rank, size = comm.getRank(), comm.getSize() + if rank == 0: + xs = np.linspace(0.0, 1.0, N + 1) + coords = np.array([[x, y] for y in xs for x in xs]) + cells = [] + for j in range(N): + for i in range(N): + a = j * (N + 1) + i + cells += [[a, a + 1, a + N + 2], [a, a + N + 2, a + N + 1]] + cells = np.asarray(cells, dtype=PETSc.IntType) + else: + coords = np.zeros((0, 2)) + cells = np.zeros((0, 3), dtype=PETSc.IntType) + dm = PETSc.DMPlex().createFromCellList(2, cells, coords, interpolate=True, comm=comm) + dm.createLabel("Walls") + dm.markBoundaryFaces("Walls", 1) + dm.labelComplete(dm.getLabel("Walls")) + part = dm.getPartitioner() + part.setType(PETSc.Partitioner.Type.SHELL) + if rank == 0: + cS, cE = dm.getHeightStratum(0) + vS, vE = dm.getDepthStratum(0) + + def walls(c): + return sum(dm.getSupportSize(e) == 1 for e in dm.getCone(c)) + + def shared_vertices(c): + return all(len(_cells_of(dm, v)) > 1 + for v in dm.getTransitiveClosure(c)[0] if vS <= v < vE) + + if where == "interior": + chosen = [c for c in range(cS, cE) if walls(c) == 0] + else: + chosen = [c for c in range(cS, cE) if walls(c) == 1 and shared_vertices(c)] + c0 = chosen[len(chosen) // 2] + rest = np.array_split(np.asarray([c for c in range(cS, cE) if c != c0]), size - 1) + part.setShellPartition(size, [1] + [len(r) for r in rest], + [c0] + [int(c) for r in rest for c in r]) + else: + part.setShellPartition(size, [0] * size, []) + dm.distribute() + return dm + + +@pytest.fixture(scope="module", params=["interior", "boundary"]) +def mesh(request): + mesh = uw.discretisation.Mesh(_owned_nothing_dm(request.param), degree=1, qdegree=2) + mesh._case = request.param + return mesh + + +def _all(flag): + return uw.mpi.comm.allreduce(bool(flag), op=MPI.LAND) + + +def _cell_vertex_coords(mesh): + """Each local cell's vertex coordinates, read from the coordinate section.""" + dm = mesh.dm + vS, vE = dm.getDepthStratum(0) + sec = dm.getCoordinateSection() + x = dm.getCoordinatesLocal().array + cS, cE = dm.getHeightStratum(0) + out = [] + for c in range(cS, cE): + pts = dm.getTransitiveClosure(c)[0] + out.append(np.array([x[sec.getOffset(v):sec.getOffset(v) + 2] + for v in pts if vS <= v < vE])) + return out + + +def _sorted_rows(a): + a = np.round(np.asarray(a), 12) + return a[np.lexsort(a.T[::-1])] + + +def _rank_zero_owns_nothing(mesh, where="interior"): + """Rank 0 holds one cell and owns no vertex, nor (interior cell) any edge, as + one verdict for every rank: a rank-0-only assert would leave the others + waiting in the next collective until the timeout.""" + cS, cE = mesh.dm.getHeightStratum(0) + leaves = set(np.asarray(mesh.dm.getPointSF().getGraph()[1]).tolist()) + owned = [] + for d in (0, 1): + pS, pE = mesh.dm.getDepthStratum(d) + owned.append(sum(1 for p in range(pS, pE) if p not in leaves)) + expected = [0, 0] if where == "interior" else [0, 1] + return _all(uw.mpi.rank != 0 or (cE - cS == 1 and owned == expected)) + + +def test_rank_zero_holds_a_cell_and_owns_nothing(mesh): + """The case under test.""" + assert _rank_zero_owns_nothing(mesh, mesh._case) + + +def _expected_nodes(verts, degree, t=None): + """A triangle's Lagrange nodes: its vertices, then (degree 2) its edge + midpoints, or (degree 3) two points on each edge, at ``t`` and ``1 - t`` + from either end, and the centroid.""" + a, b = verts[[0, 1, 2]], verts[[1, 2, 0]] + if degree == 1: + return verts + if degree == 2: + return np.vstack([verts, 0.5 * (a + b)]) + return np.vstack([verts, (1 - t) * a + t * b, t * a + (1 - t) * b, + verts.mean(axis=0, keepdims=True)]) + + +@pytest.mark.parametrize("degree", [1, 2, 3]) +def test_node_coordinates_on_every_rank(mesh, degree): + """Each cell's nodes are its vertices and, by degree, its edge midpoints, + or two points per edge and its centroid. Degree 3 places its edge points by + PETSc's node rule (``node_endpoints`` False), so ``t`` is read from the + first cell and must then hold on every edge of every cell, on every rank.""" + coords = mesh._get_coords_for_basis(degree, True) + rows = mesh._cell_node_indices(degree, True) + cells = _cell_vertex_coords(mesh) + t = None + if degree == 3 and cells: + v, nodes = cells[0], coords[rows[0]] + on_edge = [np.dot(x - v[0], v[1] - v[0]) / np.dot(v[1] - v[0], v[1] - v[0]) + for x in nodes if abs(np.cross(v[1] - v[0], x - v[0])) < 1e-12] + t = min(s for s in on_edge if s > 1e-12) + t_all = uw.mpi.comm.allgather(t) + ok = all(abs(s - t_all[0]) < 1e-12 for s in t_all if s is not None) if degree == 3 else True + for k, verts in enumerate(cells): + got = coords[rows[k]] + ok &= np.allclose(_sorted_rows(got), _sorted_rows(_expected_nodes(verts, degree, t)), + atol=1e-12) + assert _all(ok) + + +def test_projection_between_degrees_on_every_rank(mesh): + """P1 to P2 of a linear field and P2 to P1 of a quadratic one are exact at the + nodes.""" + s1 = uw.discretisation.MeshVariable("s1_1107", mesh, 1, degree=1) + s2 = uw.discretisation.MeshVariable("s2_1107", mesh, 1, degree=2) + x1 = np.asarray(s1.coords) + x2 = np.asarray(s2.coords) + s1.array[...] = (x1[:, 0] + 2.0 * x1[:, 1]).reshape(-1, 1, 1) + s2.array[...] = (x2[:, 0] ** 2).reshape(-1, 1, 1) + + up = project_to_degree(s1, target_degree=2) + down = project_to_degree(s2, target_degree=1) + n2 = mesh._get_coords_for_basis(2, True) + n1 = mesh._get_coords_for_basis(1, True) + ok = (np.allclose(up[:, 0], n2[:, 0] + 2.0 * n2[:, 1], atol=1e-12) + and np.allclose(down[:, 0], n1[:, 0] ** 2, atol=1e-12)) + assert _all(ok) + + +@pytest.mark.parametrize("where", ["interior", "boundary"]) +def test_a_degree_two_mesh_builds_on_every_rank(where): + """``Mesh(degree=2)`` gives the plex a degree-2 coordinate space, carried + across from the vertex coordinates: the vertices stay on the grid and every + edge's coordinate is its midpoint.""" + mesh = uw.discretisation.Mesh(_owned_nothing_dm(where), degree=2, qdegree=3) + assert _rank_zero_owns_nothing(mesh, where) # the case under test + dm = mesh.dm + sec = dm.getCoordinateSection() + x = dm.getCoordinatesLocal().array + + def at(p): + return x[sec.getOffset(p):sec.getOffset(p) + 2] + + ok = True + vS, vE = dm.getDepthStratum(0) + for v in range(vS, vE): + ok &= sec.getDof(v) == 2 and np.allclose(at(v) * N, np.round(at(v) * N), atol=1e-12) + eS, eE = dm.getDepthStratum(1) + for e in range(eS, eE): + a, b = dm.getCone(e) + ok &= sec.getDof(e) == 2 and np.allclose(at(e), 0.5 * (at(a) + at(b)), atol=1e-12) + assert _all(ok) + + +def _poisson_on(dm, name): + """A degree-2 Poisson problem, zero on the walls with a unit source.""" + mesh = uw.discretisation.Mesh(dm.clone(), degree=1, qdegree=3, boundaries=_WALLS) + T = uw.discretisation.MeshVariable(name, mesh, 1, degree=2) + p = uw.systems.Poisson(mesh, T) + p.constitutive_model = uw.constitutive_models.DiffusionModel + p.constitutive_model.Parameters.diffusivity = 1.0 + p.f = 1.0 + p.add_dirichlet_bc(0.0, "Walls") + p.petsc_options["ksp_rtol"] = 1.0e-10 + return p + + +def test_custom_fmg_over_a_level_where_a_rank_owns_nothing(): + """Custom-P FMG reads every level's node coordinates. Its coarse level is the + grid on which rank 0 owns nothing, under the grid's refinement. The build + installs multigrid, the solve agrees with GAMG's, and it takes the + iterations it takes on a partition where every rank owns DOFs. The + transfer's input, each level's node coordinates, is pinned exactly by the + tests above.""" + base = _owned_nothing_dm() + fine = base.refine() + coarse = uw.discretisation.Mesh(base.clone(), degree=1, qdegree=3, boundaries=_WALLS) + assert _rank_zero_owns_nothing(coarse) # the case under test + + reference = _poisson_on(fine, "Tg_1107") + reference.preconditioner = "gamg" + reference.solve() + + solver = _poisson_on(fine, "Tm_1107") + custom_mg.set_custom_fmg(solver, [coarse], builder="barycentric") + solver.solve() + + assert solver.snes.getKSP().getPC().getType() == "mg" + assert solver.snes.getConvergedReason() > 0 + assert "custom_mg.build" not in solver.pc_fallbacks + assert solver.snes.getKSP().getIterationNumber() == FMG_ITERATIONS + u, r = np.asarray(solver.Unknowns.u.array), np.asarray(reference.Unknowns.u.array) + num = uw.mpi.comm.allreduce(float(np.sum((u - r) ** 2))) + den = uw.mpi.comm.allreduce(float(np.sum(r ** 2))) + assert (num / den) ** 0.5 < 1.0e-6 diff --git a/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py new file mode 100644 index 000000000..7930a39e4 --- /dev/null +++ b/tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py @@ -0,0 +1,84 @@ +"""A transfer builder that fails on one rank is retried on every rank (#886). + +``build_transfers`` retries a failed barycentric build with the RBF builder, and +degrades to the solver's default preconditioner when that fails too. The verdict +was each rank's own: a rank whose builder raised retried alone while the others +went on into the transfer's assembly, and the ranks met in different collectives +(on kaiju at np 32 and 52, after #902's assembly failure, the run hung). + +The barycentric builder is rank-local work, so its verdict can be agreed right +after it: a failure on any rank is raised on every rank, and the retry is taken +together. Here the builder raises on rank 1 only. Contract: every rank records +the RBF rescue, geometric multigrid is installed, and the solve converges. A +check on allreduced sizes fails on every rank at once; it degrades every rank +to the default preconditioner together, as in serial. Before the fix the +ranks hung in the transfer's assembly. The thread timeout stops a rank waiting +in Python; a rank held inside a PETSc collective keeps the GIL, so a hang +there is ended by the MPI supervisor's silence limit instead. + + mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1108_custom_mg_builder_failure_collective_parallel.py +""" +import pytest + +import underworld3 as uw +from underworld3.utilities import custom_mg + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, pytest.mark.tier_b, + pytest.mark.timeout(300, method="thread")] + + +def _fmg_poisson(name): + """A single-field Poisson asking for geometric multigrid over the mesh's own + refinement hierarchy (the opportunistic, degradable build).""" + box = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, refinement=2, + regular=True, qdegree=3) + T = uw.discretisation.MeshVariable(name, box, 1, degree=2) + p = uw.systems.Poisson(box, T) + p.constitutive_model = uw.constitutive_models.DiffusionModel + p.constitutive_model.Parameters.diffusivity = 1.0 + p.add_dirichlet_bc(0.0, "Bottom") + p.add_dirichlet_bc(1.0, "Top") + p.preconditioner = "fmg" + return p + + +def test_a_builder_failure_on_one_rank_is_retried_on_every_rank(monkeypatch): + barycentric = custom_mg._BUILDERS["barycentric"] + + def fails_on_rank_one(coarse_coords, fine_coords): + if uw.mpi.rank == 1: + raise RuntimeError("synthetic barycentric failure on rank 1 (test probe)") + return barycentric(coarse_coords, fine_coords, repair_orphans=False) + + monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", fails_on_rank_one) + p = _fmg_poisson("T1108") + with pytest.warns(UserWarning, match="retrying with the 'rbf' builder"): + p.solve() + + assert p.snes.getConvergedReason() > 0 + assert p.snes.getKSP().getPC().getType() == "mg" + rec = p.pc_fallbacks["custom_mg.transfer_builder"] + assert rec["reason"] == "build_failed" and "rank 1" in rec["detail"], rec + assert "custom_mg.build" not in p.pc_fallbacks + + +def test_a_failure_every_rank_raises_degrades_every_rank(monkeypatch): + """The finest map's size is checked against the operator on allreduced sizes, + so every rank raises together. Shifting each rank's owned range by one makes + the check fail everywhere: the build is retried, then degraded to the + default preconditioner on every rank, and the solve converges on it.""" + check = custom_mg.CustomMGHierarchy._assert_finest_matches_operator + + def off_by_one(solver, finest_map, parallel): + return check(solver, finest_map._replace(rend=finest_map.rend + 1), parallel) + + monkeypatch.setattr(custom_mg.CustomMGHierarchy, "_assert_finest_matches_operator", + staticmethod(off_by_one)) + p = _fmg_poisson("T1108b") + with pytest.warns(UserWarning, match="custom-P FMG build failed"): + p.solve() + + assert p.snes.getConvergedReason() > 0 + assert p.snes.getKSP().getPC().getType() == "gamg" + assert p.pc_fallbacks["custom_mg.build"]["reason"] == "build_failed"