# Remove constant pressure solution using nullspace with PETSc

When solving (Navier-)Stokes equations with BC (Dirichlet or Neumann) on the velocity only, the pressure is only defined up to an additive constant.

There are examples in the FF documentation that use a penalization term to compensate for this issue, but I am interested in a more exact approach. I am specifically trying to use the `nullspace` argument in a PETSc Mat object to remove the constant part of the pressure.

``````real[int,int] nullsp(J.n,1);
Wh def(pconst) = [0,0,1]; // Define the nullspace (constant pressure field)
ChangeNumbering(J, pconst[], nullsp(:,0)); // convert to PETSc numbering
set(J, sparams = "-pc_type lu", nullspace = nullsp); // attach nullspace to Mat
``````

However, when Iâ€™ve tried to implement this in a MWE (attached), I do not get expected behavior. Does anyone have any insights into why this isnâ€™t working?
navier-stokes-2d-PETSc-cavity.edp (3.8 KB)
*This example is derived from `examples/hpddm/navier-stokes-2d-PETSc.edp`.

Is your RHS in the image of `J`?

Not exactly. `J` is built from a linearization of the RHS. However, since the pressure appears only as a linear term, that part of the RHS is in the image of `J`.

One issue is that you are using an exact LU factorization to deal with `J`, but since the problem is singular, the factorization may be inaccurate. You usually want to switch to field split to deal with this more efficiently from a numerical point of view.

True. I hadnâ€™t thought of that issue. So does this mean that `nullspace` only works for iterative methods? I guess that addressing the issue by projecting the problem onto the non-singular subspace (outside of the defined nullspace) would destroy the sparsity structure of the operatorâ€¦

It sounds like an easier fix would be to replace one row of the problem with a Dirichlet constraint on the pressure. However, I am not sure how to do this in parallel. Here is a MWE that works for 1 proc, but fails with >1 proc because the part of the row stored on other procs is not set to zero. Do you have any idea about how to correctly implement this? Even for 1 proc, this implementation is probably far from optimalâ€¦
navier-stokes-2d-PETSc-MWE.edp (4.2 KB)

You could multiply `pconst[]` by `J.D` to make sure that only a value on the interior of the subdomain is locked (and thus, the row will only be set by the current process).

1 Like

This helps, but I think it still doesnâ€™t communicate that the row should be set to zero everywhere except on the diagonal. I only get quadratic convergence with 1 proc.
navier-stokes-2d-PETSc-MWE2.edp (4.2 KB)

OK, letâ€™s backup a bit. What do you want to do, precisely? If you are just looking at imposing a mean pressure of 0, without penalizing the pressure block, you could simply add a Lagrange multiplier defined as

``````varf vL([u1, u2, p], [v1, v2, q])  = int2d(Th)(q);
``````

Then, your Jacobian would become

``````Mat N = [[J,  L],
[L', 0]];
``````

Doesnâ€™t this work?

Yes, you are right. Thatâ€™s exactly what I have normally done up until now, and it does work.

This detour came about because I was curious if the same result could be achieved without the extra step of augmenting the system with a Lagrange multiplier/constraint. But I am now convinced that simply locking a single DOF (common in the finite-difference literature) does not work, since the finite element method leverages a weak formulation. In a finite element approach, significant, spatially-localized errors appear near the replaced DOF since the weak statement cannot enforce both continuity and a pressure constraint simultaneously. This should have been obvious â€“ sorry to bother!

Right, locking is often done in computational solid mechanics, but Iâ€™ve never seen applied to CFD. Iâ€™m sorry itâ€™s one of these times where the solution of a difficult question is not simple and easy.

1 Like

Is it any better if you just replace p with dp/dx and introduce dp/dy?
Now you are solving for two variables but there may be benefits and how
many more non-zero matrix elements do you get?

Interesting idea, but I donâ€™t think that would be simpler or cleaner in my case.