diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 414393ed..3d521845 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/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 diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 75a42ea9..96b82241 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