Insert essential boundary values in the rotated solve's field copy-back

src/underworld3/utilities/rotated_bc.py

@@ -386,23 +386,19 @@ def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge)
             q.destroy()
             removed = True
 
-    # scatter U → velocity/pressure fields
-    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)
-
-    # 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
-    # downstream consumers (var.data / var.array / checkpoint / stats) don't read a
-    # stale value; and mark the mesh local vector stale.
-    solver.mesh._stale_lvec = True
-    for name, var in solver.fields.items():
-        target_var = getattr(var, "_base_var", var)
-        if hasattr(target_var, "_sync_lvec_to_gvec"):
-            target_var._sync_lvec_to_gvec()
-        if hasattr(target_var, "_canonical_data"):
-            target_var._canonical_data = None
+    # scatter U → velocity/pressure fields, completing each field's essential
+    # (Dirichlet) DOFs. Those are absent from the global vector, so a plain
+    # per-field scatter leaves them at ZERO — silently wrong wherever the datum
+    # g != 0 (an inhomogeneous Dirichlet wall next to a rotated boundary, #497).
+    # The insertion must run on the FULL dm, where the auxiliary vector lives:
+    # a sub-DM insertion segfaults when a BC datum references another
+    # MeshVariable. _scatter_global_to_fields does exactly that (boundary FEM on
+    # the parent local vector, then the per-field split) and carries the same
+    # cache-invalidation tail as the native post-solve copy-back.
+    # TODO(BUG): DMPlexSNESComputeBoundaryFEM inserts at time=PETSC_MIN_REAL, so
+    # a mesh.t-dependent essential datum is written as garbage — same defect on
+    # the native copy-back path. See issue #410.
+    solver._scatter_global_to_fields(U)
     return removed
 
 

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