# Remove constant pressure solution using nullspace with PETSc

**URL:** <https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712>\
**Category:** General Discussion\
**Created:** [September 18, 2023, 10:31am UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712 "2023-09-18T10:31:39Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [September 18, 2023, 10:31am UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/1 "2023-09-18T10:31:39Z")

</div>

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.

```auto
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](https://community.freefem.org/uploads/short-url/hqJvUApQPaNnhMTDnLKAn04apky.edp) (3.8 KB)  
\*This example is derived from `examples/hpddm/navier-stokes-2d-PETSc.edp`.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 18, 2023, 10:51am UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/2 "2023-09-18T10:51:20Z")

</div>

Is your RHS in the image of `J`?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [September 18, 2023, 11:15am UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/3 "2023-09-18T11:15:55Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 18, 2023, 11:38am UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/4 "2023-09-18T11:38:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [September 18, 2023, 1:31pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/5 "2023-09-18T13:31:00Z")

</div>

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](https://community.freefem.org/uploads/short-url/7OHW8wTqYRLVJazYWVjxC9BmMfq.edp) (4.2 KB)

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 18, 2023, 2:12pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/6 "2023-09-18T14:12:41Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [September 18, 2023, 2:40pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/7 "2023-09-18T14:40:36Z")

</div>

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](https://community.freefem.org/uploads/short-url/k1oNvhYTXfrh9kGhDrgo6Dc540E.edp) (4.2 KB)

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 18, 2023, 4:10pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/8 "2023-09-18T16:10:25Z")

</div>

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

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

```

Then, your Jacobian would become

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

```

Doesn’t this work?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [September 18, 2023, 4:48pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/9 "2023-09-18T16:48:43Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 18, 2023, 4:54pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/10 "2023-09-18T16:54:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 18, 2023, 6:57pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/11 "2023-09-18T18:57:46Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [September 18, 2023, 8:06pm UTC](https://community.freefem.org/t/remove-constant-pressure-solution-using-nullspace-with-petsc/2712/12 "2023-09-18T20:06:10Z")

</div>

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