Stokes PETSC solver but in a doubled-size space

Ah, I see. Thanks for clarifying. Unfortunately, I do not know about the well-posedness. You are enforcing some type of Robin-like condition resulting from the sum of the implicit Neumann conditions from integrating the viscous and pressure terms by parts, plus the vect(vh)'*vect(v) term. I’m not sure how this behaves.