Inconsistent velocity L2 norm across CPU counts in boundary-integral branch

[gthyagi]

Hi @lmoresi, I am still unable to reproduce consistent L2 norm values for velocity and pressure with the recent boundary-integral branch.

Here is the benchmark script I am using:
https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/benchmarks/annulus/Ex_Stokes_Kramer_latest.py

Could you please check if there might be any trivial mistakes in the setup and whether the boundary integrals are being called correctly?

L2 norms:

[gthyagi]

What is the issue (in plain English)?

The solver is applying a force on an internal interface (a line inside the mesh).
When the mesh is partitioned across CPUs, that interface can sometimes be represented slightly differently.

On some CPU counts (e.g., 4, 6, 8), one small interface segment is treated as one-sided instead of two-sided, meaning its contribution is either missed or handled differently.

This small missing contribution slightly alters the final solution, which causes the L2 error to vary depending on the number of CPUs.

Why the L2 norm differs across CPUs:
Parallel partitioning changes how one internal segment is assembled, so the computation is not fully partition-invariant yet.

---

What tests were done to conclude this?

Several controlled tests were performed:

1. Reproducer runs

Script: internal_bc_rank_reproducer.py
Runs were executed with different MPI sizes (-np).

Results:

• Low-error group: np = 1, 7, 10
  • v_l2_int ≈ 2.323e-03
• High-error group: np = 4, 6, 8
  • v_l2_int ≈ 6.3e-03

---

2. Label consistency checks

Compared labels:

• Internal
• UW_Boundaries

Result: labels match, so this is not a label mismatch issue.

---

3. Internal facet topology diagnostics

Examined support size per facet.

• Low-error runs:

  • 806 proper two-sided internal facets

• High-error runs:

  • 805 two-sided facets
  • 1 problematic one-sided facet copy

---

4. Side-selection A/B test

Tested internal interface side selection:

• part = 0
• part = 1

Results:

• part = 1 produced very large errors (~1.96e-01) even at np=1.
• This confirms the correct physics is on part = 0.
• The issue cannot be solved simply by flipping the side.

---

5. Fallback and overlap experiments

Tested fallback and overlap-style assembly for the problematic one-sided facet.

Result: the issue persisted.

---

Best plan of action to debug/fix next

1. Add assembly-level logging

Add logging in the Underworld core internal BC assembly to track:

• which facets are used,
• which facets are skipped,
• reasons for skipping (e.g., support != 2, ghost facet, part mismatch).

---

2. Compare working vs failing runs

Compare logs for:

• np = 7 (correct result)
• np = 8 (incorrect result)

Goal: isolate the exact facet assembly path that differs.

---

3. Patch internal facet assembly

Make internal-facet assembly ownership-consistent for distributed meshes:

• Normal two-sided facets: use the current assembly path.
• One-sided distributed artifacts: apply a deterministic fallback path so the contribution is assembled once and consistently.

---

4. Re-run validation matrix

Re-run solver with:
np = 1, 4, 6, 7, 8, 10

Requirement: L2 norms should be nearly identical across all runs.

---

5. Add a regression test

Introduce a regression test so that this issue cannot silently reappear in future changes.

[gthyagi]

I’m still unable to find a fix for the bug. Codex seems to be going in circles and hasn’t been able to debug it. I might need help from another human brain — and maybe a different AI tool 😅.

This appears to be what the recent PETSc patch is intended to address.

[gthyagi]

Short Answer

Partly, but not fully.

I reviewed plexfem-ghost-facet-fix.patch and PETSC_MR_DESCRIPTION.md.

---

What the Patch Does

The patch filters ghost facets (SF leaves) in PETSc boundary residual/integral paths to avoid double-counting shared facets.

This behavior is related to the issue we are seeing.

---

Additional Nuance Observed

Current test results show an extra complication:

• On problematic runs (np = 8), one internal facet appears with support = 1 (one-sided) on both the owner and the ghost.
• I emulated the ghost filtering behavior (keeping only the one-sided owner and skipping the ghost).
• However, the error remained high:
  v_l2_int ≈ 6.277e-03
  rather than matching the low-error branch.

---

Implication

This suggests that ghost filtering alone is not sufficient for this case.

In addition to ghost filtering, we likely need correct assembly handling for one-sided internal facets.

---

Conclusion

The PETSc patch is moving in the right direction, but it likely does not completely resolve the rank-dependent behavior currently observed.

[gthyagi]

Fixed the bug!

[gthyagi]

check #78 to issue this issue.

[lmoresi]

Fixed by PR #78 (PETSc internal-boundary ownership patch), merged via PR #70.