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.
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).
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
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.
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?
That’s exactly what I have normally done up until now, and it does work.