From 72b292d1087776938ac9c3707c5cd2a4539d2f3b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 23 Jun 2026 15:34:57 +1000 Subject: [PATCH] Nitsche BC: scale penalty by a local per-cell mesh size MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Nitsche free-slip/constraint penalty term gamma*mu/h previously used a single GLOBAL scalar h = mesh.get_min_radius() (the smallest cell anywhere). The viscosity mu was already local; only h was wrong. On a non-uniform, adaptive, or deforming mesh this mis-scales the boundary penalty: the global minimum cell size is applied on every facet (over-stiffening coarser boundary cells), and it drifts as refinement or distortion changes the global minimum. Replace it with a LOCAL, per-cell size via a new deformation-tracking field Mesh.cell_size(): - a cell-constant (degree-0, discontinuous) scalar MeshVariable holding each cell's characteristic length (self._radii), so the boundary kernel sees the adjacent cell's size; - refreshed on mesh.deform() and mesh.adapt() (and via a reinit callback in the remesh transaction), mirroring the boundary_normal machinery, so it is never stale on a moved/re-refined mesh; - the fill is purely rank-local (no collective), so it is parallel-safe. Wired into add_nitsche_bc on both SNES_Stokes_SaddlePt and SNES_Vector behind local_h=True (default). local_h=False restores the exact legacy global-h path. On a uniform mesh local == global, so existing tuned behaviour is preserved. Validation: - Existing Nitsche suite green serial + np2 (test_1060/0056/1011/1014/1061; the only np2 failure, test_1061, segfaults on unmodified development too — its LU/mumps block-constrained path uses no Nitsche). - New test_1065_nitsche_local_h (serial + np2): the penalty size is local (~6x the global minimum at a coarse free-slip boundary), tracks deformation, and still solves free-slip correctly. - Held-lid adaptive free-surface convection (fs_convection_goal4, res 24, rho_g 7.5e5, gamma 10): global-h spikes vhmax to 2873 (2 glitches, spurious h_inf mountain 2.7e-2); local-h holds vhmax at 342 (0 glitches, h_inf 7.5e-3). Medians identical => same cost, spikes removed. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 60 +++- .../discretisation/discretisation_mesh.py | 103 +++++++ tests/test_1065_nitsche_local_h.py | 282 ++++++++++++++++++ 3 files changed, 429 insertions(+), 16 deletions(-) create mode 100644 tests/test_1065_nitsche_local_h.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 1db6079e..830d6950 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -2887,7 +2887,7 @@ class SNES_Vector(SolverBaseClass): self.petsc_options["ksp_atol"] = self._tolerance * 1.0e-6 - def add_nitsche_bc(self, boundary, g=None, direction=None, gamma=10.0, theta=1): + def add_nitsche_bc(self, boundary, g=None, direction=None, gamma=10.0, theta=1, local_h=True): r"""Add Nitsche weak enforcement of a velocity constraint along a direction. For vector solvers (no pressure field), this constrains @@ -2906,6 +2906,11 @@ class SNES_Vector(SolverBaseClass): Dimensionless stabilisation parameter. theta : {-1, 0, 1}, default=1 Symmetry parameter (1=symmetric, -1=skew-symmetric). + local_h : bool, default=True + Scale the penalty by a local per-cell mesh size + (:meth:`Mesh.cell_size`) rather than the global minimum + (:meth:`Mesh.get_min_radius`). See + ``SNES_Stokes_SaddlePt.add_nitsche_bc`` for details. Warnings -------- @@ -2952,12 +2957,18 @@ class SNES_Vector(SolverBaseClass): g = sympy.Integer(0) constraint = u_dot_d - g - # Mesh size - h = uw.function.expression( - r"h_{\mathrm{Nitsche}}", - mesh.get_min_radius(), - "Nitsche mesh size parameter", - ) + # Mesh size for the penalty term (gamma*mu/h). Default: a LOCAL, + # per-cell size (mesh.cell_size()) that tracks deformation/adaptation + # so the stabilisation is scaled correctly on a non-uniform mesh. Set + # local_h=False for the legacy single global-minimum scalar. + if local_h: + h_sym = mesh.cell_size() + else: + h_sym = uw.function.expression( + r"h_{\mathrm{Nitsche}}", + mesh.get_min_radius(), + "Nitsche mesh size parameter (global)", + ).sym # Viscosity from constitutive model mu = self.constitutive_model.viscosity @@ -2974,7 +2985,7 @@ class SNES_Vector(SolverBaseClass): # f0_bd: velocity boundary residual (value term) f0_components = [] for c in range(dim): - f0_c = (gamma * mu / h.sym) * constraint * d[c] # penalty + f0_c = (gamma * mu / h_sym) * constraint * d[c] # penalty f0_c -= t_d * d[c] # consistency f0_components.append(f0_c) @@ -4689,7 +4700,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # BC = namedtuple('EssentialBC', ['components', 'fn', 'boundary', 'boundary_label_val', 'type', 'PETScID']) # self.essential_p_bcs.append(BC(components, sympy_fn, boundary, -1, 'essential', -1)) - def add_nitsche_bc(self, boundary, g=None, direction=None, normal=None, gamma=10.0, theta=1, mask=None): + def add_nitsche_bc(self, boundary, g=None, direction=None, normal=None, gamma=10.0, theta=1, mask=None, local_h=True): r"""Add Nitsche weak enforcement of a velocity constraint along a direction. Nitsche's method provides a variationally consistent alternative to @@ -4738,6 +4749,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): boundaries. Use a DG MeshVariable that is 1 on the active side and 0 on the inactive side. The mask multiplies all Nitsche terms so that only the active-side cell contributes. + local_h : bool, default=True + Scale the penalty term :math:`\gamma\mu/h` by a **local**, + per-cell mesh size (:meth:`Mesh.cell_size`, deformation- and + adaptation-tracking) rather than the single **global** minimum + cell size (:meth:`Mesh.get_min_radius`). On a non-uniform or + adaptive mesh the local size scales the stabilisation correctly + on every facet; on a uniform mesh the two coincide. Set ``False`` + to restore the legacy global-h behaviour exactly. Examples -------- @@ -4815,12 +4834,21 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # n_dot_d is always positive when n = d (it's |n|²), # so sign of PETSc face normal doesn't affect this term - # Mesh size (global estimate via UWexpression constant) - h = uw.function.expression( - r"h_{\mathrm{Nitsche}}", - mesh.get_min_radius(), - "Nitsche mesh size parameter", - ) + # Mesh size for the penalty term (gamma*mu/h). Default: a LOCAL, + # per-cell characteristic size (mesh.cell_size()), so the Nitsche + # stabilisation is correctly scaled on every facet of a non-uniform + # or adaptively-refined mesh — the boundary kernel sees the adjacent + # cell's size. The field tracks mesh deformation/adaptation. Set + # local_h=False to restore the legacy single global-minimum scalar + # (mesh.get_min_radius()); on a uniform mesh the two coincide. + if local_h: + h_sym = mesh.cell_size() + else: + h_sym = uw.function.expression( + r"h_{\mathrm{Nitsche}}", + mesh.get_min_radius(), + "Nitsche mesh size parameter (global)", + ).sym # Viscosity from constitutive model mu = self.constitutive_model.viscosity @@ -4839,7 +4867,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # = penalty + consistency + pressure flux f0_components = [] for c in range(dim): - f0_c = (gamma * mu / h.sym) * constraint * d[c] # penalty + f0_c = (gamma * mu / h_sym) * constraint * d[c] # penalty f0_c -= t_d * d[c] # consistency f0_c += p_sym * n_dot_d * d[c] # pressure flux f0_components.append(f0_c) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 32b90128..7628e503 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2595,6 +2595,91 @@ def _assemble_boundary_normal(self, var, name): accum[nonzero] /= mag[nonzero, numpy.newaxis] var.data[...] = accum + def cell_size(self): + """Local, per-cell characteristic mesh size as a scalar field symbol. + + Returns the ``.sym`` of a cell-constant (degree-0, discontinuous) + scalar MeshVariable holding each cell's characteristic length (the + ``volume**(1/dim)`` equivalent radius, i.e. ``self._radii``). Unlike + the single *global* scalar from :meth:`get_min_radius` (the smallest + cell anywhere), this varies cell to cell, so a stabilisation that + scales as :math:`1/h` — e.g. the Nitsche free-slip penalty + :math:`\\gamma\\mu/h` — is correctly scaled on every facet of a + non-uniform or adaptively-refined mesh rather than using the global + minimum (which over-penalises coarse cells and drifts as refinement + changes the global min). + + On a boundary integral the kernel sees the value of the cell adjacent + to the facet. The field is cached and rebuilt lazily; its data is + refreshed when the mesh deforms or is adapted (see :meth:`deform`), + so it tracks a moving / re-refined mesh — a stale size on a deformed + mesh would re-introduce the mis-scaling. + + On a uniform mesh every cell is the same size, so this reduces to the + global ``get_min_radius`` value everywhere and existing behaviour is + preserved to tolerance. + + Returns + ------- + sympy scalar + The cell-size field symbol, for use in JIT-compiled residuals. + """ + var = self._cell_size_var() + return var.sym[0] + + def _cell_size_var(self): + """Lazily create / fetch the per-cell size MeshVariable (filled). + + Mirrors the per-boundary normal machinery (:meth:`boundary_normal`): + a small ``reinit`` MeshVariable owned by the mesh, refreshed from the + current geometry. The reinit callback re-fills it during any remesh + transaction (deform / mover sweep), and :meth:`deform` / :meth:`adapt` + re-fill it explicitly so BCs that captured ``cell_size()`` at setup + read the new geometry at solve time. + """ + import underworld3 as uw + + if getattr(self, "_cell_size_variable", None) is None: + existing = self.vars.get("_h_cell") + var = existing if existing is not None else uw.discretisation.MeshVariable( + "_h_cell", self, 1, degree=0, continuous=False, + remesh_policy="reinit") + self._cell_size_variable = var + + def _refresh(): + try: + self._assemble_cell_size(var) + except Exception: + pass + + var._remesh_reinit_callback = _refresh + + self._assemble_cell_size(self._cell_size_variable) + return self._cell_size_variable + + def _assemble_cell_size(self, var): + """Fill ``var`` (degree-0 scalar) with each cell's characteristic size. + + Uses the per-cell characteristic lengths ``self._radii`` computed by + :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0 + discontinuous variable's local DOFs and ``self._radii`` are BOTH + indexed by this rank's cell-stratum order, so a direct assignment is + correct on every rank. + + This is deliberately a purely RANK-LOCAL operation (no ``var.coords`` + access, no collective): mixing a rank-local fast path with a + collective fallback would diverge across ranks and deadlock, because + ``var.coords`` triggers the collective ``_get_coords_for_basis``.""" + radii = numpy.asarray(self._radii).reshape(-1) + # Empty partition (no local cells): nothing to fill on this rank. + if radii.size == 0 or var.data.shape[0] == 0: + return + # Assign over the common length. In practice these match exactly (same + # local cell set / ordering); the slice only guards a stray off-by-ghost + # mismatch without ever taking a collective path on a subset of ranks. + n = min(var.data.shape[0], radii.shape[0]) + var.data[:n, 0] = radii[:n] + @property def Gamma_P1(self): """Projected P1 boundary normals as a sympy Matrix. @@ -3051,6 +3136,13 @@ def _do_move(): self._assemble_boundary_normal(_var, _nm) except Exception: pass + # Likewise refresh the local cell-size field (Nitsche penalty scaling) + # so its cell-constant data tracks the deformed geometry. + if getattr(self, "_cell_size_variable", None) is not None: + try: + self._assemble_cell_size(self._cell_size_variable) + except Exception: + pass return result def _deform_mesh(self, new_coords: numpy.ndarray, verbose=False, @@ -6057,6 +6149,17 @@ def mesh_update_callback(array, change_context): # Note: Surfaces were already notified at the start of adapt() # They will lazily recompute distance fields when accessed + # Refresh the local cell-size field from the adapted geometry. The + # generic variable transfer above re-interpolates it (meaningless for + # a geometric size); re-fill it from the new mesh's per-cell radii so + # the Nitsche penalty scales correctly after re-refinement. + if getattr(self, "_cell_size_variable", None) is not None: + try: + self._assemble_cell_size(self._cell_size_variable) + except Exception as e: + if verbose: + print(f"[{uw.mpi.rank}] Warning: cell-size refresh failed: {e}", flush=True) + # Mark solvers for rebuild for solver in self._equation_systems_register: if solver is not None and hasattr(solver, '_rebuild_after_mesh_update'): diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py new file mode 100644 index 00000000..31484294 --- /dev/null +++ b/tests/test_1065_nitsche_local_h.py @@ -0,0 +1,282 @@ +"""Local-h scaling for the Nitsche free-slip penalty. + +``add_nitsche_bc`` stabilises the constraint with a term ``gamma*mu/h``. +Historically ``h`` was a single GLOBAL scalar (``mesh.get_min_radius()`` — the +smallest cell anywhere). On a non-uniform / adaptive / deforming mesh that +mis-scales the boundary penalty: the *global minimum* cell size is applied on +every facet, even where the boundary cells are much coarser, and it drifts as +refinement (or distortion) changes the global minimum. ``local_h=True`` (the +default) replaces it with a per-cell, deformation-tracking size field +(``mesh.cell_size()``). + +What these tests assert — and why they are framed around the *mechanism* +rather than a boundary-``v.n`` leak comparison: + + On a STATIC graded mesh the global-h penalty is the global *minimum*, so it + is never weaker than the correct local penalty — it *over*-stiffens the + coarse boundary cells. Over-stiffening gives a SMALLER ``v.n`` leak there + (it approaches a hard constraint), so a naive "local enforces v.n better + than global" leak test is actually backwards on a static mesh. The real + costs of global-h are (a) the wrong asymptotic conditioning of the velocity + block on graded meshes and (b) spurious boundary-velocity spikes on a + *deforming* mesh once a distorted/refined cell drags the global minimum far + below the boundary cell size. (b) is exercised by the held-lid free-surface + integration in ``~/+Simulations/fs_convection_goal4`` (vhmax ~3000 -> ~300). + +So here we assert the things that are both TRUE and decisive on a static mesh: + 1. the penalty size is genuinely LOCAL (matches each cell's size, not the + global minimum) — and at a coarse free-slip boundary it is many times the + global minimum that global-h would wrongly use; + 2. the field TRACKS deformation (it is not stale after ``mesh.deform``) — + this is the property whose absence caused the free-surface spike; + 3. local-h still solves free-slip correctly (matches the essential-BC + reference and keeps ``v.n`` small on both the coarse and fine ends). + +Run: pixi run python -m pytest tests/test_1065_nitsche_local_h.py -v +""" + +import os +from enum import Enum + +import numpy as np +import pytest +import sympy +from mpi4py import MPI + +import underworld3 as uw +from underworld3.discretisation import Mesh +from underworld3.coordinates import CoordinateSystemType + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +# Parallel-safe global reductions. Every assertion below is made on a value +# that is IDENTICAL on all ranks (a global reduction), so the ranks never +# diverge on a pass/fail decision — a rank-local assertion that passes on one +# rank and fails on another would leave the passing rank spinning forever in +# the next collective. +def _gmin(a): + a = np.asarray(a).reshape(-1) + loc = float(a.min()) if a.size else float("inf") + return uw.mpi.comm.allreduce(loc, op=MPI.MIN) + + +def _gmax(a): + a = np.asarray(a).reshape(-1) + loc = float(a.max()) if a.size else float("-inf") + return uw.mpi.comm.allreduce(loc, op=MPI.MAX) + + +def _any(flag): + return bool(uw.mpi.comm.allreduce(bool(flag), op=MPI.LOR)) + + +class _boundaries_2D(Enum): + Bottom = 11 + Top = 12 + Right = 13 + Left = 14 + + +class _boundary_normals_2D(Enum): + Bottom = sympy.Matrix([0, 1]) + Top = sympy.Matrix([0, 1]) + Right = sympy.Matrix([1, 0]) + Left = sympy.Matrix([1, 0]) + + +def _graded_box(h_fine=0.04, h_coarse=0.12): + """Unit box graded FINE at the Bottom edge, COARSE at the Top edge. + + The Top free-slip boundary therefore sits on cells ~``h_coarse`` while the + global minimum cell size is ~``h_fine`` (down at the Bottom), so global-h + would over-stiffen the Top penalty by ~``h_coarse/h_fine``. + """ + os.makedirs(".meshes", exist_ok=True) + fname = f".meshes/uw_test_graded_box_{h_fine}_{h_coarse}.msh" + + if uw.mpi.rank == 0: + import gmsh + + gmsh.initialize() + gmsh.option.setNumber("General.Verbosity", 0) + gmsh.model.add("GradedBox") + p1 = gmsh.model.geo.add_point(0.0, 0.0, 0.0, meshSize=h_fine) + p2 = gmsh.model.geo.add_point(1.0, 0.0, 0.0, meshSize=h_fine) + p3 = gmsh.model.geo.add_point(0.0, 1.0, 0.0, meshSize=h_coarse) + p4 = gmsh.model.geo.add_point(1.0, 1.0, 0.0, meshSize=h_coarse) + l1 = gmsh.model.geo.add_line(p1, p2, tag=_boundaries_2D.Bottom.value) + l2 = gmsh.model.geo.add_line(p2, p4, tag=_boundaries_2D.Right.value) + l3 = gmsh.model.geo.add_line(p4, p3, tag=_boundaries_2D.Top.value) + l4 = gmsh.model.geo.add_line(p3, p1, tag=_boundaries_2D.Left.value) + cl = gmsh.model.geo.add_curve_loop((l1, l2, l3, l4)) + surface = gmsh.model.geo.add_plane_surface([cl]) + gmsh.model.geo.synchronize() + for l, b in ((l1, "Bottom"), (l2, "Right"), (l3, "Top"), (l4, "Left")): + gmsh.model.add_physical_group(1, [l], l) + gmsh.model.set_physical_name(1, l, b) + gmsh.model.addPhysicalGroup(2, [surface], 99999) + gmsh.model.setPhysicalName(2, 99999, "Elements") + gmsh.model.mesh.generate(2) + gmsh.write(fname) + gmsh.finalize() + + # ensure the file is on disk before non-root ranks read it (cold cache) + uw.mpi.barrier() + + return Mesh( + fname, + degree=1, + qdegree=3, + boundaries=_boundaries_2D, + boundary_normals=_boundary_normals_2D, + coordinate_system_type=CoordinateSystemType.CARTESIAN, + useMultipleTags=True, + useRegions=True, + markVertices=True, + ) + + +def _box_wobble(X0, amp): + """Smooth interior perturbation tapering to zero on all four box edges, + so boundaries stay put and the mesh stays valid for small ``amp``.""" + c = np.asarray(X0).copy() + bump = np.sin(np.pi * c[:, 0]) * np.sin(np.pi * c[:, 1]) + c[:, 0] += amp * bump + c[:, 1] += amp * bump + return c + + +# -------------------------------------------------------------------------- +# 1. The penalty size is LOCAL, not the global minimum +# -------------------------------------------------------------------------- +def test_cell_size_is_local_per_cell(): + """``mesh.cell_size()`` is a per-cell field equal to each cell's + characteristic size (``mesh._radii``), not the single global minimum.""" + mesh = _graded_box() + h = mesh.cell_size() # sympy symbol -> backed by a P0 field + field = np.asarray(mesh._cell_size_variable.data[:, 0]).reshape(-1) + radii = np.asarray(mesh._radii).reshape(-1) + + # field exactly mirrors the per-cell characteristic size (rank-local check, + # reduced to a single global pass/fail so all ranks agree) + n = min(field.shape[0], radii.shape[0]) + assert not _any(not np.allclose(field[:n], radii[:n])) + + # the mesh is genuinely graded (so local != global is meaningful) — use + # GLOBAL min/max so the full fine-to-coarse range is seen in parallel too + gfmin, gfmax = _gmin(field), _gmax(field) + assert gfmax / gfmin > 3.0 + + # the global scalar that global-h would use is just the minimum cell size + assert np.isclose(mesh.get_min_radius(), gfmin, rtol=1e-6) + + +def test_local_h_at_coarse_freeslip_boundary_exceeds_global_min(): + """At the COARSE Top free-slip boundary the LOCAL penalty size is many + times the global minimum. global-h would over-stiffen that penalty by + exactly this factor; local-h scales it correctly.""" + mesh = _graded_box(h_fine=0.04, h_coarse=0.12) + # build/exercise the field; its data equals mesh._radii (asserted in + # test_cell_size_is_local_per_cell), so we read the per-cell sizes directly + # from _radii / _centroids — a rank-local lookup, avoiding the collective + # arbitrary-point uw.function.evaluate (which deadlocks in parallel). + _ = mesh.cell_size() + cen = np.asarray(mesh._centroids) + radii = np.asarray(mesh._radii).reshape(-1) + + near_top = cen[:, 1] > 0.85 # cells adjacent to the Top free-slip edge + h_top = radii[near_top] + + g = mesh.get_min_radius() # global minimum (the value global-h uses) + h_top_min = _gmin(h_top) # smallest LOCAL size along the Top (global) + # The Top free-slip cells are clearly coarser than the global minimum, so + # global-h would over-stiffen the Top penalty (gamma*mu/g) by ~ h_top/g. + assert h_top_min > 2.0 * g + + +# -------------------------------------------------------------------------- +# 2. The field tracks deformation (it is not stale) — the free-surface bug +# -------------------------------------------------------------------------- +def test_cell_size_tracks_deformation(): + """After ``mesh.deform`` the cell-size field is refreshed to the new + geometry — a stale size on a deformed mesh is exactly what re-introduced + the Nitsche mis-scaling on the free surface.""" + mesh = _graded_box() + _ = mesh.cell_size() + h_before = mesh._cell_size_variable.data[:, 0].copy() + + X = np.asarray(mesh.X.coords).copy() + moved = mesh.deform(_box_wobble(X, amp=0.04)) + assert moved # geometry actually changed + + h_after = mesh._cell_size_variable.data[:, 0].copy() + radii_after = np.asarray(mesh._radii).reshape(-1) + + # not stale: the field changed with the geometry SOMEWHERE (global OR) ... + nb = min(h_after.shape[0], h_before.shape[0]) + assert _any(not np.allclose(h_after[:nb], h_before[:nb])) + # ... and on EVERY rank it equals the freshly recomputed per-cell sizes. + na = min(h_after.shape[0], radii_after.shape[0]) + assert not _any(not np.allclose(h_after[:na], radii_after[:na])) + + +# -------------------------------------------------------------------------- +# 3. local-h still solves free-slip correctly (back-compat / correctness) +# -------------------------------------------------------------------------- +def _solve_freeslip(mesh, method, gamma=10.0): + v = uw.discretisation.MeshVariable( + "U", mesh, mesh.dim, degree=2, vtype=uw.VarType.VECTOR) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.saddle_preconditioner = 1.0 + x, y = mesh.X + stokes.bodyforce = sympy.Matrix([0, sympy.cos(sympy.pi * x)]) + stokes.add_dirichlet_bc((0.0, 0.0), "Left") + stokes.add_dirichlet_bc((0.0, 0.0), "Right") + if method == "essential": + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Top") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Bottom") + else: + local = method == "nitsche_local" + stokes.add_nitsche_bc("Top", gamma=gamma, local_h=local) + stokes.add_nitsche_bc("Bottom", gamma=gamma, local_h=local) + stokes.tolerance = 1e-8 + stokes.petsc_options["snes_type"] = "ksponly" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.solve() + return v + + +def _vn_rms(mesh, v, boundary): + Gamma = mesh.Gamma + vn = v.sym.dot(Gamma) + num = float(uw.maths.BdIntegral(mesh, fn=vn ** 2, boundary=boundary).evaluate()) + length = float(uw.maths.BdIntegral(mesh, fn=1.0, boundary=boundary).evaluate()) + return (num / length) ** 0.5 + + +def test_local_h_freeslip_solution_is_correct(): + """On the graded mesh, local-h free-slip Nitsche reproduces the + essential-BC reference and keeps v.n small on BOTH the coarse (Top) and + fine (Bottom) ends — the constraint holds everywhere, independent of the + local refinement.""" + mesh_ref = _graded_box() + v_ref = _solve_freeslip(mesh_ref, "essential") + vrms_ref = float(np.sqrt(uw.maths.Integral( + mesh_ref, v_ref.sym.dot(v_ref.sym)).evaluate())) + + mesh = _graded_box() + v = _solve_freeslip(mesh, "nitsche_local") + vrms = float(np.sqrt(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate())) + + # matches the exact free-slip (hard-constraint) solution + assert abs(vrms - vrms_ref) / vrms_ref < 0.02 + + # constraint enforced on both the coarse and the fine boundary + top = _vn_rms(mesh, v, "Top") + bot = _vn_rms(mesh, v, "Bottom") + assert top < 1.0e-3 + assert bot < 1.0e-3