Branch Review — knepley/feature-friction-assembly

Contents

1. Scope of the diff
2. CRITICAL — local vector passed to DMProjectFunction()
3. CRITICAL — hybrid key part hard-coded to 0/1/2
4. HIGH — nondeterministic side assignment for cohesive keys
5. MEDIUM — no test coverage for the new/changed code paths
6. Verified correct / not findings
7. Suppressed LOW findings
8. Suggested fixes (diffs)

Scope of the diff

Findings

PetscCall(DMGetLocalVector(dm, &locb));
PetscCall(VecSet(locb, PETSC_NINFINITY));
PetscCall(DMProjectFunction(dm, 0., lfuncs, lctxs, INSERT_VALUES, locb));   // ← locb is local
PetscCall(DMPlexInsertBounds(dm, PETSC_TRUE, 0., locb));

keys[0].part  = 0;
keys[1].part  = 1;
keys[2].part  = 2;

PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);

for (PetscInt i = 0; i < bdf0Nk; ++i) {
  if (bdf0keys[i].field == 0) {
    if (!label0) { label0 = ...; value0A = bdf0keys[i].value; value0B = value0A; }
    else if (bdf0keys[i].value != value0A) value0B = bdf0keys[i].value;
  }
  ...
}

Suggested fixes (diffs)

Fix B — preserve registered part, make cohesive side ordering deterministic

Addresses CRITICAL 2 and HIGH.

--- a/src/snes/utils/dmplexsnes.c
+++ b/src/snes/utils/dmplexsnes.c
@@ -420,6 +420,7 @@ PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec locX, Vec locF, PetscCtx
       IS       cohesiveCells;
       DMLabel  label0 = NULL, label1 = NULL;
       PetscInt value0A = 0, value0B = 0, value1 = 0;
+      PetscInt part0A = 0, part0B = 0, part1 = 0;

       // In the future, we need a way to construct the keys for both sides and the surface, and also the cellIS from the label
       for (PetscInt i = 0; i < bdf0Nk; ++i) {
@@ -427,28 +428,41 @@ PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec locX, Vec locF, PetscCtx
           if (!label0) {
             label0  = bdf0keys[i].label;
             value0A = bdf0keys[i].value;
+            part0A  = bdf0keys[i].part;
             value0B = value0A;
+            part0B  = part0A;
           } else if (bdf0keys[i].value != value0A) {
             value0B = bdf0keys[i].value;
+            part0B  = bdf0keys[i].part;
           }
         }
         if (bdf0keys[i].field == 1) {
           label1 = bdf0keys[i].label;
           value1 = bdf0keys[i].value;
+          part1  = bdf0keys[i].part;
         }
       }
+      // Order sides by value so the key assignment does not depend on PetscWeakFormGetKeys() hash map iteration order
+      if (value0B < value0A) {
+        PetscInt tmpValue = value0A, tmpPart = part0A;
+
+        value0A = value0B;
+        part0A  = part0B;
+        value0B = tmpValue;
+        part0B  = tmpPart;
+      }
       keys[0].label = label0;
       keys[0].value = value0A;
       keys[0].field = 0;
-      keys[0].part  = 0;
+      keys[0].part  = part0A;
       keys[1].label = label0;
       keys[1].value = value0B;
       keys[1].field = 0;
-      keys[1].part  = 1;
+      keys[1].part  = part0B;
       keys[2].label = label1;
       keys[2].value = value1;
       keys[2].field = 1;
-      keys[2].part  = 2;
+      keys[2].part  = part1;
       PetscCall(CreateSurfaceCellIS_Private(plex, &cohesiveCells));
       PetscCall(DMPlexComputeResidualHybridByKey(plex, keys, cohesiveCells, PETSC_MIN_REAL, locX, NULL, 0.0, locF, ctx));
       PetscCall(ISDestroy(&cohesiveCells));

part now reflects whatever value was actually used when the residual functions were registered (normally 0), so PetscWeakFormGetBdResidual()'s exact-match lookup succeeds for all three keys instead of only keys[0]. Sorting by value makes the side assignment reproducible instead of depending on PetscHMapFormGetKeys()'s pointer-hash-derived iteration order — and it happens to match the convention already used in ex69.c (negative side = lower label value, positive side = higher). This does not address the unused-bdf1keys gap noted above, which needs a design decision rather than a mechanical fix.

Fix A — use DMProjectFunctionLocal() instead of DMProjectFunction()

Addresses CRITICAL 1.

--- a/src/snes/utils/dmplexsnes.c
+++ b/src/snes/utils/dmplexsnes.c
@@ -1229,7 +1243,7 @@ PetscErrorCode DMPlexSetSNESVariableBounds(DM dm, SNES snes)

     PetscCall(DMGetLocalVector(dm, &locb));
     PetscCall(VecSet(locb, PETSC_NINFINITY));
-    PetscCall(DMProjectFunction(dm, 0., lfuncs, lctxs, INSERT_VALUES, locb));
+    PetscCall(DMProjectFunctionLocal(dm, 0., lfuncs, lctxs, INSERT_VALUES, locb));
     PetscCall(DMPlexInsertBounds(dm, PETSC_TRUE, 0., locb));
     PetscCall(DMLocalToGlobalBegin(dm, locb, INSERT_VALUES, lb));
     PetscCall(DMLocalToGlobalEnd(dm, locb, INSERT_VALUES, lb));
@@ -1242,7 +1256,7 @@ PetscErrorCode DMPlexSetSNESVariableBounds(DM dm, SNES snes)

     PetscCall(DMGetLocalVector(dm, &locb));
     PetscCall(VecSet(locb, PETSC_INFINITY));
-    PetscCall(DMProjectFunction(dm, 0., ufuncs, uctxs, INSERT_VALUES, locb));
+    PetscCall(DMProjectFunctionLocal(dm, 0., ufuncs, uctxs, INSERT_VALUES, locb));
     PetscCall(DMPlexInsertBounds(dm, PETSC_FALSE, 0., locb));
     PetscCall(DMLocalToGlobalBegin(dm, locb, INSERT_VALUES, ub));
     PetscCall(DMLocalToGlobalEnd(dm, locb, INSERT_VALUES, ub));

Both call sites just swap the function name; DMProjectFunction() and DMProjectFunctionLocal() share the same argument list, differing only in whether the final Vec is global or local.