From d82b7f8cd60f34081b46420bbbfd8bfc03e09474 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 4 Aug 2026 18:49:08 +1000 Subject: [PATCH] Insert essential boundary values in the rotated solve's field copy-back Essential-BC DOFs are absent from the global vector, so the rotated path's per-field globalToLocal scatter left them at ZERO in the output fields wherever the Dirichlet datum was non-zero. The solve itself was always right -- the residual path inserts the values internally -- but every field-based diagnostic (projection, Integral, evaluate, renders) read a garbage boundary strip. Measured: far-field sigma_xy of 0.79 for a true 1.00, and a projected stress row off by 3x one fault length out. Invisible to every earlier rotated test because their Dirichlet walls were homogeneous: zero was accidentally the datum. Fix: petsc_dm_insert_boundary_values in cython/petsc_discretisation.pyx -- wrapping DMPlexInsertBoundaryValues, which petsc4py does not expose; the same insertion the consistent-boundary-flux paths already use (issues #407/#411) -- called per field after the scatter in _finalize_rotated_solution. The regression test drives an inhomogeneous lid with rotated free-slip sides and asserts the FIELD carries the wall values. The tell that found it: the divergence theorem, Integral(2 e_xy) against its own boundary term. The rotated boundary's prescribed datum and the sigma_nn reaction recovery were never affected (both live in the global vector). Fixes #497. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9 --- .../cython/petsc_discretisation.pyx | 17 ++++++++ src/underworld3/utilities/rotated_bc.py | 10 ++++- tests/test_1018_rotated_freeslip.py | 41 +++++++++++++++++++ 3 files changed, 67 insertions(+), 1 deletion(-) diff --git a/src/underworld3/cython/petsc_discretisation.pyx b/src/underworld3/cython/petsc_discretisation.pyx index cd5a248f0..a15613057 100644 --- a/src/underworld3/cython/petsc_discretisation.pyx +++ b/src/underworld3/cython/petsc_discretisation.pyx @@ -309,3 +309,20 @@ def petsc_vec_concatenate( inputVecs ): def petsc_get_swarm_coord_name( sdm ): return + + +def petsc_dm_insert_boundary_values(dm_, lvec_, time=0.0): + """Insert essential boundary values into a LOCAL vector of ``dm_``. + + Constrained DOFs are absent from the global system, so a global-to-local + scatter alone leaves them at zero wherever the essential datum g != 0 — + the same gap the consistent-boundary-flux paths close internally (issues + #407/#411). A solve that bypasses SNES's own copy-back (the rotated + strong-BC path) calls this to complete its output fields. Wraps + ``DMPlexInsertBoundaryValues``, which petsc4py does not expose. + """ + cdef DM dm = dm_ + cdef Vec lv = lvec_ + cdef PetscReal t = time + CHKERRQ(DMPlexInsertBoundaryValues(dm.dm, PETSC_TRUE, lv.vec, t, + NULL, NULL, NULL)) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 414393ed1..efa38d94b 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -386,11 +386,19 @@ def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) q.destroy() removed = True - # scatter U → velocity/pressure fields + # scatter U → velocity/pressure fields. Constrained (essential-BC) DOFs + # are absent from the global vector, so the scatter leaves them at ZERO in + # the local field — silently wrong wherever the datum g != 0 (an + # inhomogeneous Dirichlet wall next to a rotated boundary). Complete each + # field with the DS's own essential values, exactly as the native SNES + # copy-back does (the #407/#411 insertion, via the cython shim). + from underworld3.cython.petsc_discretisation import \ + petsc_dm_insert_boundary_values for name, var in solver.fields.items(): sg = U.getSubVector(solver._subdict[name][0]) solver._subdict[name][1].globalToLocal(sg, var.vec) U.restoreSubVector(solver._subdict[name][0], sg) + petsc_dm_insert_boundary_values(solver._subdict[name][1], var.vec) # Parity with the normal solve's post-scatter sync (pyx: after the field copy-back): # refresh the enhanced-variable gvec cache and drop the canonical-data cache so diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 75a42ea92..96b82241f 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -726,3 +726,44 @@ def test_rotated_freeslip_nonlinear_prescribed_normal_datum(): err = np.abs(vn - target).max() assert err < 1e-8, f"nonlinear u.n=cos(theta) not imposed: max nodal error {err:.2e}" assert vn.max() > 0.9 and vn.min() < -0.9, "prescribed normal velocity magnitude wrong" + + +def test_rotated_solve_fields_carry_inhomogeneous_dirichlet_walls(): + """The copy-back gap: essential DOFs are absent from the global vector, so + the rotated path's field scatter left them at ZERO wherever the Dirichlet + datum g != 0 — the solve was right, every field-based diagnostic + (projection, integral, evaluate) read a garbage boundary strip. Caught by + the split-fault work (far-field stress off by 20%); fixed by the + DMPlexInsertBoundaryValues shim in the copy-back. Homogeneous walls hid + this from every earlier rotated test — zero happens to be their datum. + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + x, y = mesh.X + v = uw.discretisation.MeshVariable("vIB", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("pIB", mesh, 1, degree=1, + continuous=False) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + s.tolerance = 1e-8 + s.petsc_use_pressure_nullspace = True + # Inhomogeneous Dirichlet lid and floor, rotated free-slip sides: the + # combination that exposes the gap. + s.add_dirichlet_bc((y - 0.5, 0.0), "Top") + s.add_dirichlet_bc((y - 0.5, 0.0), "Bottom") + s.add_rotated_freeslip_bc(0, "Left") + s.add_rotated_freeslip_bc(0, "Right") + s.solve() + + vc = np.asarray(v.coords) + vd = np.asarray(v.data) + for name, mask, target in ( + ("Top", vc[:, 1] > 1 - 1e-9, +0.5), + ("Bottom", vc[:, 1] < 1e-9, -0.5)): + assert mask.sum() > 0 + err = np.abs(vd[mask, 0] - target).max() + assert err < 1e-10, ( + f"{name} wall u_x in the FIELD is off by {err:.2e}; the rotated " + "copy-back dropped the inhomogeneous essential values") + assert np.abs(vd[mask, 1]).max() < 1e-10