Discontinuities at MPI partition boundaries with manual periodic MPC on large high-order curved mesh

dolfinx
,
3

Thank you @dokken for clarifying that the explicit ghostUpdate + backsubstitution after solve() is redundant since it’s already done internally.

I added those calls during diagnostics following advice in thread 12401 where they seemed needed in older versions. Removing them does not change the symptom — the partition-boundary discontinuities persist.

Setup recap (since this is the core question):

  • Poisson problem ∇²u = f on a curved hex27 mesh from CFD data (~17M DOFs, P4 Lagrange)
  • Two periodic translations (pitch + span), no Dirichlet boundary except a single pin
  • ~200k MPC slaves
  • 16 MPI ranks
  • Periodicity built via manual mpc.add_constraint(slaves, masters, coeffs, owners, offsets) because the slave/master pairs come from external CFD mesh data

What I tried after your reply:

I attempted to switch from manual add_constraint to the higher-level create_periodic_constraint_topological and create_periodic_constraint_geometrical, but both immediately crash with:

RuntimeError: Newton method failed to converge for non-affine geometry

This matches the issue from thread 13511 and thread 13421, where you suggested raising the pull-back Newton tolerance (PR #104). That tolerance does not seem exposed via the Python API in v0.10.0, so I’d need to patch and rebuild the C++ side — which is non-trivial in our conda-managed cluster environment.

The manual add_constraint path is therefore the only route that works for me (no crash). With that path, the solve converges, and all internal consistency checks pass on every rank in parallel:

  • u.function_space is mpc.function_space → True
  • max_i | u[slave_i] - sum_j coeff_ij * u[master_ij] | = 0.000e+00
  • Ghost-owned-value consistency check (compare every ghost dof’s value to its owned counterpart on the owning rank, via KDTree on allgathered coords) → 0.000e+00 on every rank

Yet the solution still shows visible jumps along MPI partition boundaries when visualized in ParaView — confirmed to coincide with the partition seams by overlaying a DG-0 rank field, and the number of jumps scales with the rank count (2 jumps with n=2, 15 jumps with n=16).

My questions:

  1. Is there a way to expose the pull-back Newton tolerance via the Python API, so I could try create_periodic_constraint_topological on my curved mesh?
  2. Given that all the standard consistency checks pass on the manual add_constraint path, is there something subtle in the matrix assembly with that low-level API that could lead to partition-dependent results — even when the slave-master values are perfectly enforced?