diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index 2fe0c649e..74075e043 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/src/underworld3/cython/petsc_compat.h @@ -45,6 +45,39 @@ PetscErrorCode DMSetAuxiliaryVec_UW(DM dm, DMLabel label, PetscInt value, PetscI return DMSetAuxiliaryVec(dm, label, value, part, aux); } +PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS ds, PetscInt bd, PetscWeakForm *wf) +{ + PetscCall(PetscDSGetBoundary(ds, bd, wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); + return PETSC_SUCCESS; +} + +PetscErrorCode UW_DMPlexComputeResidualByKeyVolumeOnly(DM dm, PetscFormKey key, IS cellIS, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *ctx) +{ + DM vdm = NULL; + PetscSection section = NULL; + PetscSection global_section = NULL; + PetscDS ds = NULL; + PetscDS vds = NULL; + Vec aux = NULL; + + PetscCall(DMClone(dm, &vdm)); + PetscCall(DMGetLocalSection(dm, §ion)); + PetscCall(DMSetLocalSection(vdm, section)); + PetscCall(DMGetGlobalSection(dm, &global_section)); + PetscCall(DMSetGlobalSection(vdm, global_section)); + PetscCall(DMCopyFields(dm, PETSC_DETERMINE, PETSC_DETERMINE, vdm)); + PetscCall(DMCreateDS(vdm)); + PetscCall(DMGetDS(dm, &ds)); + PetscCall(DMGetDS(vdm, &vds)); + PetscCall(PetscDSCopyConstants(ds, vds)); + PetscCall(PetscDSCopyEquations(ds, vds)); + PetscCall(DMGetAuxiliaryVec(dm, key.label, key.value, key.part, &aux)); + if (aux) PetscCall(DMSetAuxiliaryVec(vdm, key.label, key.value, key.part, aux)); + PetscCall(DMPlexComputeResidualByKey(vdm, key, cellIS, time, locX, locX_t, t, locF, ctx)); + PetscCall(DMDestroy(&vdm)); + return PETSC_SUCCESS; +} + // copy paste function signitures from $PETSC_DIR/include/petscds.h - would be nice to automate this. #define UW_SIG_F0 PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[] #define UW_SIG_G0 PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[] diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index 23785c99c..f864e0d8c 100644 --- a/src/underworld3/cython/petsc_extras.pxi +++ b/src/underworld3/cython/petsc_extras.pxi @@ -31,10 +31,20 @@ cdef CHKERRQ(PetscErrorCode ierr): cdef int interr = ierr if ierr != 0: raise RuntimeError(f"PETSc error code '{interr}' was encountered.\nhttps://www.mcs.anl.gov/petsc/petsc-current/include/petscerror.h.html") +cdef extern from "petscdstypes.h": + ctypedef void *PetscWeakForm "PetscWeakForm" + ctypedef struct PetscFormKey: + PetscDMLabel label + PetscInt value + PetscInt field + PetscInt part + cdef extern from "petsc_compat.h": PetscErrorCode PetscDSAddBoundary_UW( PetscDM, DMBoundaryConditionType, const char[], const char[] , PetscInt, PetscInt, PetscInt *, void (*)(), void (*)(), PetscInt, const PetscInt *, void *) PetscErrorCode DMSetAuxiliaryVec_UW(PetscDM, PetscDMLabel, PetscInt, PetscInt, PetscVec) + PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS, PetscInt, PetscWeakForm *) + PetscErrorCode UW_DMPlexComputeResidualByKeyVolumeOnly( PetscDM, PetscFormKey, PetscIS, PetscReal, PetscVec, PetscVec, PetscReal, PetscVec, void *) # PetscErrorCode UW_PetscDSSetBdResidual(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, void*, PetscInt, void*) PetscErrorCode UW_PetscDSSetBdResidual(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*) @@ -52,6 +62,9 @@ cdef extern from "petsc_compat.h": cdef extern from "petsc.h" nogil: PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[]) PetscErrorCode DMPlexSNESComputeBoundaryFEM( PetscDM, void *, void *) + PetscErrorCode DMPlexSNESComputeResidualFEM( PetscDM, PetscVec, PetscVec, void *) + PetscErrorCode DMPlexComputeResidualByKey( PetscDM, PetscFormKey, PetscIS, PetscReal, PetscVec, PetscVec, PetscReal, PetscVec, void *) + PetscErrorCode DMPlexComputeBdResidualSingle( PetscDM, PetscWeakForm, PetscFormKey, PetscVec, PetscVec, PetscReal, PetscVec ) # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, void *, void *, void *) # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *) PetscErrorCode DMPlexComputeGeometryFVM( PetscDM dm, PetscVec *cellgeom, PetscVec *facegeom) @@ -62,6 +75,7 @@ cdef extern from "petsc.h" nogil: PetscErrorCode PetscDSSetJacobian( PetscDS, PetscInt, PetscInt, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn) PetscErrorCode PetscDSSetJacobianPreconditioner( PetscDS, PetscInt, PetscInt, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn) PetscErrorCode PetscDSSetResidual( PetscDS, PetscInt, PetscDSResidualFn, PetscDSResidualFn ) + PetscErrorCode PetscDSGetWeakForm( PetscDS, PetscWeakForm * ) PetscErrorCode PetscDSSetBdJacobian( PetscDS, PetscInt, PetscInt, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn) PetscErrorCode PetscDSSetBdJacobianPreconditioner( PetscDS, PetscInt, PetscInt, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn) @@ -83,6 +97,7 @@ cdef extern from "petsc.h" nogil: PetscErrorCode DMGetRegionDS(PetscDM dm, PetscDMLabel label, PetscIS *fields, PetscDS *ds, PetscDS *dsIn) PetscErrorCode DMGetRegionNumDS(PetscDM dm, PetscInt num, PetscDMLabel *label, PetscIS *fields, PetscDS *ds, PetscDS *dsIn) PetscErrorCode DMSetRegionNumDS(PetscDM dm, PetscInt num, PetscDMLabel label, PetscIS fields, PetscDS ds, PetscDS dsIn) + PetscErrorCode DMGetDS(PetscDM dm, PetscDS *ds) PetscErrorCode DMGetNumDS(PetscDM dm, PetscInt *num) PetscErrorCode DMGetCellDS(PetscDM dm, PetscInt point, PetscDS *ds, PetscDS *dsIn) PetscErrorCode PetscDSSetCoordinateDimension(PetscDS ds, PetscInt dim) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 0c71cf880..4848c7195 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -6682,6 +6682,270 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): if uw.mpi.rank == 0 and self.verbose: print(f"Region DS: inactive region '{label_name}' gets trivial DS", flush=True) + def compute_volume_residual_fields( + self, + time=None, + verbose=False, + cell_indices=None, + residual_field_id=None, + include_boundary_terms=False, + ): + """Return PETSc FEM residual fields in each solver field's local layout. + + This is a low-level diagnostic hook for post-processing derived + boundary quantities such as consistent-boundary-flux traction. By + default it calls PETSc's ``DMPlexSNESComputeResidualFEM`` directly. If + ``cell_indices`` is supplied, it instead calls + a UW wrapper around ``DMPlexComputeResidualByKey`` on a cloned DM with + a copied ``PetscDS`` that has no registered boundary objects, so the + selected-cell path returns volume terms only. Set + ``include_boundary_terms=True`` to call PETSc's original keyed + residual behavior, which appends registered boundary residuals. The + returned arrays are local to each rank and have the same flat layout as + the corresponding MeshVariable PETSc vector. + """ + cdef DM _time_dm_residual + cdef DM dm + cdef Vec xvec + cdef Vec fvec + cdef PetscFormKey key + cdef IS ccell_is + cdef PetscReal residual_time = 0.0 + cdef PetscReal implicit_form_time = -1.7976931348623157e308 + + self._build(verbose, False, None) + + if time is not None: + if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'): + t_nd = float(uw.non_dimensionalise(time)) + else: + t_nd = float(time) + _time_dm_residual = self.dm + UW_DMSetTime(_time_dm_residual.dm, t_nd) + residual_time = t_nd + + self.mesh.update_lvec() + self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self._update_constants() + + gvec = self.dm.getGlobalVec() + xlocal = self.dm.getLocalVec() + flocal = self.dm.getLocalVec() + gvec.setArray(0.0) + xlocal.setArray(0.0) + flocal.setArray(0.0) + + try: + for name, var in self.fields.items(): + sgvec = gvec.getSubVector(self._subdict[name][0]) + subdm = self._subdict[name][1] + subdm.localToGlobal(var.vec, sgvec) + gvec.restoreSubVector(self._subdict[name][0], sgvec) + + self.dm.globalToLocal(gvec, xlocal) + + dm = self.dm + xvec = xlocal + fvec = flocal + if cell_indices is None: + CHKERRQ(DMPlexSNESComputeResidualFEM(dm.dm, xvec.vec, fvec.vec, NULL)) + else: + if residual_field_id is None: + residual_field_id = 0 + cell_is = PETSc.IS().createGeneral( + list(cell_indices), comm=PETSc.COMM_SELF + ) + try: + ccell_is = cell_is + key.label = NULL + key.value = 0 + key.field = residual_field_id + key.part = 0 + if include_boundary_terms: + CHKERRQ(DMPlexComputeResidualByKey( + dm.dm, key, ccell_is.iset, implicit_form_time, + xvec.vec, NULL, residual_time, fvec.vec, NULL, + )) + else: + CHKERRQ(UW_DMPlexComputeResidualByKeyVolumeOnly( + dm.dm, key, ccell_is.iset, implicit_form_time, + xvec.vec, NULL, residual_time, fvec.vec, NULL, + )) + finally: + cell_is.destroy() + + local_section = self.dm.getLocalSection() + pStart, pEnd = local_section.getChart() + out = {} + + for name, var in self.fields.items(): + field_id = getattr(var, "_solver_field_id", None) + if field_id is None: + field_id = getattr(var, "field_id", None) + if field_id is None: + continue + + is_field = None + created_is_field = False + if name == "velocity" and getattr(self, "_velocity_is", None) is not None: + is_field = self._velocity_is + elif name == "pressure" and getattr(self, "_pressure_is", None) is not None: + is_field = self._pressure_is + elif getattr(self, "_multiplier_is", None) is not None and name in self._multiplier_is: + is_field = self._multiplier_is[name] + else: + indices = [] + for point in range(pStart, pEnd): + dof = local_section.getFieldDof(point, field_id) + if dof > 0: + offset = local_section.getFieldOffset(point, field_id) + for i in range(dof): + indices.append(offset + i) + + is_field = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) + created_is_field = True + + try: + subvec = flocal.getSubVector(is_field) + try: + out[name] = np.array(subvec.array, copy=True) + finally: + flocal.restoreSubVector(is_field, subvec) + finally: + if created_is_field: + is_field.destroy() + + return out + finally: + self.dm.restoreLocalVec(flocal) + self.dm.restoreLocalVec(xlocal) + self.dm.restoreGlobalVec(gvec) + + def compute_boundary_residual_fields(self, boundary, time=None, verbose=False, residual_field_id=0): + """Return the registered FEM boundary residual for one named boundary. + + This is a low-level diagnostic hook for weak-boundary-condition + debugging. It assembles PETSc's boundary residual terms registered on + ``boundary`` through ``DMPlexComputeBdResidualSingle``. For Nitsche + free slip, this includes the full registered weak boundary residual, + not only the scalar penalty term. The returned arrays are local to each + rank and have the same flat layout as the corresponding MeshVariable + PETSc vector. + """ + cdef DM _time_dm_boundary_residual + cdef DM dm + cdef Vec xvec + cdef Vec fvec + cdef PetscFormKey key + cdef PetscDS ds + cdef PetscWeakForm wf + cdef DMLabel c_label + + self._build(verbose, False, None) + + boundary_bc = None + for bc in self.natural_bcs: + if bc.boundary == boundary and bc.f_id == residual_field_id: + boundary_bc = bc + break + if boundary_bc is None: + raise ValueError( + f"No natural/Nitsche boundary residual is registered for " + f"boundary '{boundary}' and field {residual_field_id}." + ) + + if time is not None: + if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'): + t_nd = float(uw.non_dimensionalise(time)) + else: + t_nd = float(time) + _time_dm_boundary_residual = self.dm + UW_DMSetTime(_time_dm_boundary_residual.dm, t_nd) + + self.mesh.update_lvec() + self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self._update_constants() + + gvec = self.dm.getGlobalVec() + xlocal = self.dm.getLocalVec() + flocal = self.dm.getLocalVec() + gvec.setArray(0.0) + xlocal.setArray(0.0) + flocal.setArray(0.0) + + try: + for name, var in self.fields.items(): + sgvec = gvec.getSubVector(self._subdict[name][0]) + subdm = self._subdict[name][1] + subdm.localToGlobal(var.vec, sgvec) + gvec.restoreSubVector(self._subdict[name][0], sgvec) + + self.dm.globalToLocal(gvec, xlocal) + + dm = self.dm + xvec = xlocal + fvec = flocal + CHKERRQ(DMGetDS(dm.dm, &ds)) + CHKERRQ(UW_PetscDSGetBoundaryWeakForm( + ds, boundary_bc.PETScID, &wf, + )) + + c_label = self.dm.getLabel("UW_Boundaries") + key.label = c_label.dmlabel + key.value = boundary_bc.boundary_label_val + key.field = residual_field_id + key.part = 0 + CHKERRQ(DMPlexComputeBdResidualSingle( + dm.dm, wf, key, xvec.vec, NULL, 0.0, fvec.vec, + )) + + local_section = self.dm.getLocalSection() + pStart, pEnd = local_section.getChart() + out = {} + + for name, var in self.fields.items(): + field_id = getattr(var, "_solver_field_id", None) + if field_id is None: + field_id = getattr(var, "field_id", None) + if field_id is None: + continue + + is_field = None + created_is_field = False + if name == "velocity" and getattr(self, "_velocity_is", None) is not None: + is_field = self._velocity_is + elif name == "pressure" and getattr(self, "_pressure_is", None) is not None: + is_field = self._pressure_is + elif getattr(self, "_multiplier_is", None) is not None and name in self._multiplier_is: + is_field = self._multiplier_is[name] + else: + indices = [] + for point in range(pStart, pEnd): + dof = local_section.getFieldDof(point, field_id) + if dof > 0: + offset = local_section.getFieldOffset(point, field_id) + for i in range(dof): + indices.append(offset + i) + + is_field = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) + created_is_field = True + + try: + subvec = flocal.getSubVector(is_field) + try: + out[name] = np.array(subvec.array, copy=True) + finally: + flocal.restoreSubVector(is_field, subvec) + finally: + if created_is_field: + is_field.destroy() + + return out + finally: + self.dm.restoreLocalVec(flocal) + self.dm.restoreLocalVec(xlocal) + self.dm.restoreGlobalVec(gvec) + def _ensure_local_field_index_sets(self, clvec, local_section): """Build (once) and cache the LOCAL index sets that decompose a parent-DM local vector into the per-field MeshVariable storage: velocity, pressure diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index fb03e1b8d..36cbd6ffc 100644 --- a/src/underworld3/meshing/spherical.py +++ b/src/underworld3/meshing/spherical.py @@ -557,6 +557,13 @@ class boundaries(Enum): Internal = 12 Upper = 13 + # This generator builds a SINGLE shell volume [radiusInner, radiusOuter] with + # the radiusInternal sphere *embedded* as a conformal internal surface (the + # `Internal` boundary). Because it is one OCC volume, gmsh cannot emit Inner/ + # Outer region physical groups — so the Inner/Outer cell regions used by + # mesh.extract_region() are created AFTER import by classifying each cell by + # its centroid radius relative to radiusInternal (exact: the embedded surface + # is conformal, so no cell straddles it). See PR #242. class regions(Enum): Inner = 101 Outer = 102 @@ -571,105 +578,118 @@ class regions(Enum): else: uw_filename = filename - # Check if r_i is greater than 0 if radiusInner <= 0: raise ValueError("The inner radius must be greater than 0.") + if not radiusInner < radiusInternal < radiusOuter: + raise ValueError( + "SphericalShellInternalBoundary requires " + "radiusInner < radiusInternal < radiusOuter." + ) if uw.mpi.rank == 0: gmsh.initialize() gmsh.option.setNumber("General.Verbosity", gmsh_verbosity) gmsh.model.add("SphereShell_with_Internal_Surface") - # Create three concentric spheres and use OCC fragment to split - # into two non-overlapping shell volumes sharing the internal surface - ball_outer = gmsh.model.occ.addSphere(0, 0, 0, radiusOuter) - ball_internal = gmsh.model.occ.addSphere(0, 0, 0, radiusInternal) - ball_inner = gmsh.model.occ.addSphere(0, 0, 0, radiusInner) + # Create the spherical shell volume. + outer = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusOuter) + inner = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInner) + gmsh.model.occ.cut( + [(3, outer)], + [(3, inner)], + removeObject=True, + removeTool=True, + ) - # Fragment creates non-overlapping pieces from the boolean intersection - out_dimtags, out_map = gmsh.model.occ.fragment( - [(3, ball_outer)], - [(3, ball_internal), (3, ball_inner)], + # Create an internal shell only to obtain a clean spherical surface at + # radiusInternal. That surface is embedded into the shell volume below; + # the duplicate volume and duplicate lower surface are removed before + # meshing. + internal = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInternal) + inner_copy = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInner) + gmsh.model.occ.cut( + [(3, internal)], + [(3, inner_copy)], + removeObject=True, + removeTool=True, ) gmsh.model.occ.synchronize() gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cellSize) - # Identify volumes and surfaces by bounding box - # For a sphere, bbox diagonal = sqrt(3) * radius - volumes = gmsh.model.getEntities(3) - surfaces = gmsh.model.getEntities(2) - def bbox_radius(dimtag): - """Estimate the sphere radius from a bounding box diagonal.""" + """Estimate a concentric sphere radius from the bounding box.""" bb = gmsh.model.get_bounding_box(dimtag[0], dimtag[1]) return np.sqrt(bb[3]**2 + bb[4]**2 + bb[5]**2) / np.sqrt(3.0) - inner_vols = [] - outer_vols = [] - solid_ball_vols = [] # r < radiusInner — to be removed - - for vol in volumes: - r_est = bbox_radius(vol) - if np.isclose(r_est, radiusInner, atol=cellSize): - solid_ball_vols.append(vol) - elif np.isclose(r_est, radiusInternal, atol=cellSize): - inner_vols.append(vol) - elif np.isclose(r_est, radiusOuter, atol=cellSize): - outer_vols.append(vol) - - # Remove the solid inner ball (r < radiusInner) - if solid_ball_vols: - gmsh.model.occ.remove(solid_ball_vols, recursive=True) - gmsh.model.occ.synchronize() - - # Re-query after removal volumes = gmsh.model.getEntities(3) - surfaces = gmsh.model.getEntities(2) + shell_vols = [ + vol + for vol in volumes + if np.isclose(bbox_radius(vol), radiusOuter, atol=cellSize * 0.5) + ] + duplicate_vols = [ + vol + for vol in volumes + if np.isclose(bbox_radius(vol), radiusInternal, atol=cellSize * 0.5) + ] + + if len(shell_vols) != 1: + raise RuntimeError( + "Could not identify the spherical-shell volume while building " + "SphericalShellInternalBoundary." + ) + + shell_vol = shell_vols[0] + shell_boundary = { + dimtag[1] + for dimtag in gmsh.model.getBoundary([shell_vol], oriented=False, recursive=False) + if dimtag[0] == 2 + } - # Classify surfaces by bounding box radius - for surface in surfaces: + lower_surface_tags = [] + internal_surface_tags = [] + upper_surface_tags = [] + duplicate_lower_tags = [] + + for surface in gmsh.model.getEntities(2): + surface_tag = surface[1] r_est = bbox_radius(surface) if np.isclose(r_est, radiusInner, atol=cellSize * 0.5): - gmsh.model.addPhysicalGroup( - surface[0], [surface[1]], - boundaries.Lower.value, name=boundaries.Lower.name, - ) + if surface_tag in shell_boundary: + lower_surface_tags.append(surface_tag) + else: + duplicate_lower_tags.append(surface_tag) elif np.isclose(r_est, radiusOuter, atol=cellSize * 0.5): - gmsh.model.addPhysicalGroup( - surface[0], [surface[1]], - boundaries.Upper.value, name=boundaries.Upper.name, - ) + upper_surface_tags.append(surface_tag) elif np.isclose(r_est, radiusInternal, atol=cellSize * 0.5): - gmsh.model.addPhysicalGroup( - surface[0], [surface[1]], - boundaries.Internal.value, name=boundaries.Internal.name, - ) - - # Classify remaining volumes into Inner and Outer - inner_vol_tags = [v[1] for v in inner_vols if v not in solid_ball_vols] - outer_vol_tags = [v[1] for v in outer_vols] - # Re-classify from current volumes in case tags changed after removal - inner_vol_tags = [] - outer_vol_tags = [] - for vol in volumes: - r_est = bbox_radius(vol) - if r_est < radiusInternal + cellSize * 0.5: - inner_vol_tags.append(vol[1]) - else: - outer_vol_tags.append(vol[1]) - - # Region physical groups - if inner_vol_tags: - gmsh.model.addPhysicalGroup(3, inner_vol_tags, - regions.Inner.value, name=regions.Inner.name) - if outer_vol_tags: - gmsh.model.addPhysicalGroup(3, outer_vol_tags, - regions.Outer.value, name=regions.Outer.name) - - # Combined elements group - all_vol_tags = inner_vol_tags + outer_vol_tags - gmsh.model.addPhysicalGroup(3, all_vol_tags, 99999, "Elements") + internal_surface_tags.append(surface_tag) + + if not lower_surface_tags or not upper_surface_tags or not internal_surface_tags: + raise RuntimeError( + "Could not identify Lower, Internal, and Upper spherical surfaces " + "while building SphericalShellInternalBoundary." + ) + + gmsh.model.mesh.embed(2, internal_surface_tags, shell_vol[0], shell_vol[1]) + + remove_dimtags = duplicate_vols + [(2, tag) for tag in duplicate_lower_tags] + if remove_dimtags: + gmsh.model.remove_entities(remove_dimtags, recursive=False) + gmsh.model.occ.remove(remove_dimtags, recursive=False) + gmsh.model.occ.synchronize() + + gmsh.model.addPhysicalGroup( + 2, lower_surface_tags, boundaries.Lower.value, name=boundaries.Lower.name + ) + gmsh.model.addPhysicalGroup( + 2, internal_surface_tags, boundaries.Internal.value, name=boundaries.Internal.name + ) + gmsh.model.addPhysicalGroup( + 2, upper_surface_tags, boundaries.Upper.value, name=boundaries.Upper.name + ) + + gmsh.model.addPhysicalGroup(shell_vol[0], [shell_vol[1]], 99999, "Elements") gmsh.model.mesh.generate(3) gmsh.write(uw_filename) @@ -730,6 +750,28 @@ def spherical_mesh_refinement_callback(dm): new_mesh.regions = regions + # Materialise the Inner/Outer cell regions on the DM (single-volume embed + # design — see note above — so gmsh emitted no region groups). Classify each + # LOCAL cell by centroid radius vs radiusInternal: the embedded internal + # surface is conformal, so every cell lies wholly inside or outside it, and + # the test is geometric/rank-local (partition-independent — no SF needed). + # extract_region("Inner"/"Outer") filters the DM on these labels. + _region_dm = new_mesh.dm + _cStart, _cEnd = _region_dm.getHeightStratum(0) + for _rname in (regions.Inner.name, regions.Outer.name): + if _region_dm.hasLabel(_rname): + _region_dm.removeLabel(_rname) + _region_dm.createLabel(_rname) + _inner_label = _region_dm.getLabel(regions.Inner.name) + _outer_label = _region_dm.getLabel(regions.Outer.name) + for _c in range(_cStart, _cEnd): + _centroid = _region_dm.computeCellGeometryFVM(_c)[1] + _rc = float(np.linalg.norm(_centroid)) + if _rc < radiusInternal: + _inner_label.setValue(_c, regions.Inner.value) + else: + _outer_label.setValue(_c, regions.Outer.value) + # Full spherical shell with internal boundary: 3 rigid rotation modes x, y, z = new_mesh.X new_mesh._nullspace_rotations = [ diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py index 0cd9e18d9..9080fd9ad 100644 --- a/tests/test_0502_boundary_integrals.py +++ b/tests/test_0502_boundary_integrals.py @@ -293,6 +293,59 @@ def test_bd_integral_annulus_internal_normal_tangential(): assert abs(value) < 0.05, f"Expected ~0, got {value}" +# --- Spherical shell internal boundary tests --- + +from underworld3.meshing import SphericalShellInternalBoundary + +_R_SHELL_INNER = 0.55 +_R_SHELL_INTERNAL = 0.775 +_R_SHELL_OUTER = 1.0 +_mesh_spherical_internal = None + + +def _get_spherical_internal_mesh(): + global _mesh_spherical_internal + if _mesh_spherical_internal is None: + _mesh_spherical_internal = SphericalShellInternalBoundary( + radiusOuter=_R_SHELL_OUTER, + radiusInternal=_R_SHELL_INTERNAL, + radiusInner=_R_SHELL_INNER, + cellSize=0.25, + degree=1, + qdegree=2, + ) + uw.discretisation.MeshVariable( + "T_spherical_internal", _mesh_spherical_internal, 1, degree=1 + ) + return _mesh_spherical_internal + + +@pytest.mark.level_2 +@pytest.mark.tier_b +def test_bd_integral_spherical_internal_boundary_areas(): + """SphericalShellInternalBoundary preserves Lower/Internal/Upper labels. + + Overrides the module-level level_1/tier_a marks: this builds a full 3D + gmsh+embed mesh (not a seconds-scale level_1 op), and the embed generator + is not yet production-soaked for tier_a. See PR #242 review. + """ + + mesh_spherical = _get_spherical_internal_mesh() + expected_areas = { + "Lower": 4.0 * np.pi * _R_SHELL_INNER**2, + "Internal": 4.0 * np.pi * _R_SHELL_INTERNAL**2, + "Upper": 4.0 * np.pi * _R_SHELL_OUTER**2, + } + + for boundary, expected in expected_areas.items(): + value = uw.maths.BdIntegral(mesh_spherical, fn=1.0, boundary=boundary).evaluate() + relative_error = abs(value - expected) / expected + assert relative_error < 0.06, ( + f"{boundary} area should be close to {expected:.4f}; " + f"got {value:.4f} (relative error {relative_error:.3f})" + ) + + def _build_spherical_shell_for_integrals(): from underworld3.meshing import SphericalShell diff --git a/tests/test_1064_constrained_spherical_shell_response.py b/tests/test_1064_constrained_spherical_shell_response.py new file mode 100644 index 000000000..56acda5ad --- /dev/null +++ b/tests/test_1064_constrained_spherical_shell_response.py @@ -0,0 +1,234 @@ +"""3-D spherical-shell constrained free-slip response regression. + +This test records facts exposed by the Zhong et al. (2008)-style benchmark: + +* the validated Nitsche/default field-split path reproduces the Zhong velocity + scale for this low-resolution response case; +* a direct-LU diagnostic path gives matching Nitsche and constrained responses, + but it does not reproduce the validated Nitsche/default response and should + not be treated as the benchmark reference; +* the practical fast grouped-Schur constrained path currently does not reproduce + the Zhong velocity response and remains an expected failure. + +Run: + pixi run -e amr-dev pytest -q tests/test_1064_constrained_spherical_shell_response.py +""" + +from functools import cache + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_3, pytest.mark.slow, pytest.mark.tier_c] + +RADIUS_INNER = 0.55 +RADIUS_INTERNAL = 0.775 +RADIUS_OUTER = 1.0 +CELL_SIZE = 1.0 / 8.0 +HARMONIC_DEGREE = 2 +NITSCHE_GAMMA = 10.0 + +ZHONG_SURFACE_VELOCITY = 1.006e-2 +ZHONG_CMB_VELOCITY = 1.186e-2 + + +@cache +def solve_response(method, solver_mode): + mesh = uw.meshing.SphericalShellInternalBoundary( + radiusOuter=RADIUS_OUTER, + radiusInternal=RADIUS_INTERNAL, + radiusInner=RADIUS_INNER, + cellSize=CELL_SIZE, + qdegree=2, + degree=1, + ) + + velocity = uw.discretisation.MeshVariable( + f"U_{method}_{solver_mode}", + mesh, + mesh.dim, + degree=2, + vtype=uw.VarType.VECTOR, + ) + pressure = uw.discretisation.MeshVariable( + f"P_{method}_{solver_mode}", + mesh, + 1, + degree=1, + continuous=True, + ) + + theta = mesh.CoordinateSystem.xR[1] + unit_r = mesh.CoordinateSystem.unit_e_0 + y_l0 = sympy.assoc_legendre(HARMONIC_DEGREE, 0, sympy.cos(theta)) + harmonic_norm = 4.0 * np.pi / (2 * HARMONIC_DEGREE + 1) + + if method == "constrained": + stokes = uw.systems.Stokes_Constrained( + mesh, + velocityField=velocity, + pressureField=pressure, + ) + elif method == "nitsche": + stokes = uw.systems.Stokes( + mesh, + velocityField=velocity, + pressureField=pressure, + ) + else: + raise ValueError(f"Unknown response method {method!r}") + + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, 0.0, 0.0]) + stokes.add_natural_bc(y_l0 * unit_r, mesh.boundaries.Internal.name) + + if method == "nitsche": + stokes.add_nitsche_bc("Upper", normal=unit_r, gamma=NITSCHE_GAMMA) + stokes.add_nitsche_bc("Lower", normal=-unit_r, gamma=NITSCHE_GAMMA) + else: + stokes.add_constraint_bc( + "Upper", + g=0.0, + normal=unit_r, + augmentation_base=1.0e4, + degree=2, + ) + stokes.add_constraint_bc( + "Lower", + g=0.0, + normal=-unit_r, + augmentation_base=1.0e4, + degree=2, + ) + + stokes.petsc_use_nullspace = True + stokes.tolerance = 1.0e-7 + stokes.petsc_options["snes_type"] = "ksponly" + + if solver_mode == "monolithic": + stokes.petsc_options["ksp_type"] = "preonly" + stokes.petsc_options["pc_type"] = "lu" + stokes.petsc_options["pc_factor_mat_solver_type"] = "mumps" + stokes.petsc_options["pc_use_amat"] = None + elif solver_mode == "default": + pass + elif solver_mode == "fieldsplit_exact": + if method != "constrained": + raise ValueError( + "fieldsplit_exact mode is only defined for constrained runs" + ) + stokes.petsc_options["pc_fieldsplit_schur_precondition"] = "selfp" + stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "preonly" + stokes.petsc_options["fieldsplit_velocity_pc_type"] = "lu" + stokes.petsc_options["fieldsplit_1_ksp_type"] = "preonly" + stokes.petsc_options["fieldsplit_1_pc_type"] = "lu" + elif solver_mode == "fast_schur": + if method != "constrained": + raise ValueError("fast_schur mode is only defined for constrained runs") + stokes.petsc_options["pc_fieldsplit_schur_precondition"] = "selfp" + stokes.petsc_options["fieldsplit_1_ksp_type"] = "preonly" + stokes.petsc_options["fieldsplit_1_pc_type"] = "gasm" + else: + raise ValueError(f"Unknown solver mode {solver_mode!r}") + + stokes.solve() + + horizontal_v2 = velocity.sym.dot(velocity.sym) - velocity.sym.dot(unit_r) ** 2 + + surface_velocity = np.sqrt( + uw.maths.BdIntegral(mesh, horizontal_v2, boundary="Upper").evaluate() + / ( + RADIUS_OUTER**2 + * HARMONIC_DEGREE + * (HARMONIC_DEGREE + 1) + * harmonic_norm + ) + ) + cmb_velocity = np.sqrt( + uw.maths.BdIntegral(mesh, horizontal_v2, boundary="Lower").evaluate() + / ( + RADIUS_INNER**2 + * HARMONIC_DEGREE + * (HARMONIC_DEGREE + 1) + * harmonic_norm + ) + ) + + return ( + float(surface_velocity), + float(cmb_velocity), + int(stokes.snes.getConvergedReason()), + ) + + +def test_default_nitsche_matches_zhong_velocity_response(): + surface_velocity, cmb_velocity, snes_reason = solve_response( + "nitsche", + "default", + ) + + assert snes_reason > 0 + assert abs(surface_velocity - ZHONG_SURFACE_VELOCITY) / ZHONG_SURFACE_VELOCITY < 0.05 + assert abs(cmb_velocity - ZHONG_CMB_VELOCITY) / ZHONG_CMB_VELOCITY < 0.05 + + +def test_direct_lu_diagnostic_constrained_matches_direct_lu_diagnostic_nitsche(): + nitsche_surface, nitsche_cmb, nitsche_reason = solve_response( + "nitsche", + "monolithic", + ) + constrained_surface, constrained_cmb, constrained_reason = solve_response( + "constrained", + "monolithic", + ) + + assert nitsche_reason > 0 + assert constrained_reason > 0 + assert abs(constrained_surface - nitsche_surface) / nitsche_surface < 0.01 + assert abs(constrained_cmb - nitsche_cmb) / nitsche_cmb < 0.01 + + +@pytest.mark.xfail( + reason=( + "Known constrained field-split failure: LU sub-solves in the " + "velocity | [p,h] preconditioner still do not reproduce the validated " + "Nitsche/default velocity response." + ), + strict=True, +) +def test_lu_subsolve_fieldsplit_constrained_matches_default_nitsche_response(): + nitsche_surface, nitsche_cmb, nitsche_reason = solve_response( + "nitsche", + "default", + ) + constrained_surface, constrained_cmb, constrained_reason = solve_response( + "constrained", + "fieldsplit_exact", + ) + + assert nitsche_reason > 0 + assert constrained_reason > 0 + assert abs(constrained_surface - nitsche_surface) / nitsche_surface < 0.01 + assert abs(constrained_cmb - nitsche_cmb) / nitsche_cmb < 0.01 + + +@pytest.mark.xfail( + reason=( + "Known fast grouped-Schur constrained response failure for the " + "3-D SphericalShellInternalBoundary Zhong-style load." + ), + strict=True, +) +def test_fast_schur_constrained_matches_zhong_velocity_response(): + surface_velocity, cmb_velocity, snes_reason = solve_response( + "constrained", + "fast_schur", + ) + + assert snes_reason > 0 + assert abs(surface_velocity - ZHONG_SURFACE_VELOCITY) / ZHONG_SURFACE_VELOCITY < 0.05 + assert abs(cmb_velocity - ZHONG_CMB_VELOCITY) / ZHONG_CMB_VELOCITY < 0.05