# Solver for unsteady N-S equation

**URL:** <https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351>\
**Category:** General Discussion\
**Created:** [June 26, 2024, 1:05pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351 "2024-06-26T13:05:21Z")\
**Posts on this page:** 14\
**Page:** 2

<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:** [October 17, 2024, 8:19am UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/21 "2024-10-17T08:19:49Z")

</div>

> what is the proper way to solve it efficiently?

What are you using currently? And why is not giving you appropriate performance?

> will PETSc reuse the inversion of the matrix?

Yes.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [October 17, 2024, 6:51pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/22 "2024-10-17T18:51:24Z")

</div>

Unfortunately, I got a blow up for my new code, every variable looks weird even in the first time step.  
I am afraid that my boundary condition is imposed in a wrong way. It seems that constructing rhs by  
`varf rhs(u,v)=int2d(Th)(uold*v)+on(1,u=0); real[int] b=rhs(0,Vh,tgv=-1);`  
is wrong. 😓 Maybe a mass matrix multiply uold `M*uold[]` is correct, I just thought they were the same thing.  
Actually I am not sure with the time derivative term discretized by finite element.

Anyway, I will figure out what’s going wrong as soon as possible.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [June 13, 2025, 12:17pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/23 "2025-06-13T12:17:15Z")

</div>

Hello, I am working on pressure Poisson equation, and your tutorial (Section 8, example 5, Poisson equation without Dirichlet BC) has been very helpful. In the process, I came across a couple of questions:

1. The augmented system` [[A, b], [b', 0]]` introduces one additional degree of freedom. Why is this extra DoF assigned to process 0?

2. How should I set up an appropriate solver? I tried `gamg+gmres`, but the performance is poor — it’s very slow and even diverges when the mesh size reaches n=200. I think this is due to the zero entry on the diagonal.

---

<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:** [June 13, 2025, 12:36pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/24 "2025-06-13T12:36:03Z")

</div>

1. It is the convention I chose. All centralized arrays are centralized on rank #0.
2. If you have a single constraint, you can compute the exact Schur complement (of dimension 1-by-1 with `PCFIELDSPLIT` and it will converge in a single iteration).

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [June 16, 2025, 1:39pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/25 "2025-06-16T13:39:23Z")

</div>

Hi, @prj  
Thanks for your reply. I’ve tried what you suggested, and the following `sparams` gives a converge result, but I’m not sure if it is actually the correct or recommended configuration:

```auto
set(N ,sparams = "-ksp_type fgmres -pc_type fieldsplit "
                + "-pc_fieldsplit_schur_precondition self -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_detect_saddle_point "
                + "-fieldsplit_0_pc_type gamg -fieldsplit_0_ksp_type preonly ");

```

Could you please take a quick look at my solver settings?

---

<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:** [June 16, 2025, 6:53pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/26 "2025-06-16T18:53:49Z")

</div>

If they give you satisfactory results, then good. It’s not what I recommended just prior, but that is fine.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [January 27, 2026, 4:44pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/27 "2026-01-27T16:44:15Z")

</div>

Hello @prj,

It has been a few days since our last discussion.

For a Poisson problem with pure Neumann boundary conditions (i.e. PPE), is it possible to handle the null space using `MatNullSpace`? I know a similar problem was posted by @cmd, but I am not sure how to deal with it in practice.

If yes, how should it be specified through the PETSc interface?

I tried pinning one DoF via `SetBC()`, but this introduces a Dirichlet constraint, so the solution is spatially affected. I also tried removing the constant mode by modifying the RHS (subtracting the mean via an inner product), but this does not seem to work correctly in parallel.

Best regards,  
Haocheng Weng

---

<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:** [January 27, 2026, 8:41pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/28 "2026-01-27T20:41:24Z")

</div>

There is a named parameter in the `set()` command to specify the nullspace.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [January 28, 2026, 8:12am UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/29 "2026-01-28T08:12:26Z")

</div>

Hello @prj,

Thanks for the tip about the named parameter in `set()` for specifying the null space.

I searched for examples and tried the following minimal test for a Poisson problem with pure Neumann boundary conditions:

```ff++
load "PETSc"
macro dimension()2//
include "macro_ddm.idp"

macro grad(u)[dx(u), dy(u)]//
func Pk = P2;

mesh Th = square(100,100);
Mat A;
MatCreate(Th, A, Pk);

fespace Vh(Th, Pk);
func f = 1 + x - y;

varf vA(u,v) = intN(Th)( grad(u)'*grad(v) ) + intN(Th)( f*v );
matrix Loc = vA(Vh, Vh, tgv = -2);
A = Loc;

real[int] b = vA(0, Vh);

// constant nullspace
real[int,int] Rb(Vh.ndof, 1);
Rb(:,0) = 1.0;

set(A,
    sparams = "-ksp_type preonly -pc_type lu -ksp_monitor_true_residual",
    nullspace = Rb);

Vh u;
u[] = A^-1 * b;
plotD(Th, u, cmm = "Solution");

```

Solvers such as `gmres + gamg` diverge.  
Using `lu + preonly` does return a solution, but the result is distorted: the global trend looks reasonable, while local structures appear incorrect.

 ![snapshot](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/8/83c7e7b4721ac80153bac2885cf2e466a4de22a7.png)

Is specifying the null space alone sufficient in this case?  
Best regards,  
Haocheng Weng

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [January 28, 2026, 9:31am UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/30 "2026-01-28T09:31:00Z")

</div>

Dear Haocheng Weng,  
theoretically, imposing a condition u(x0)=0 gives a unique solution,  
where x0 is an arbitrary point.  
Did you try with setBC, imposing a single dof to zero? This needs to be on a single point of the whole domain, hence for example only on proc 0.

Another difficulty is that for the existence of a solution, it is necessary that \int f=0 (to see that take v=1 as test function). In particular you need to subtract to f its average on the whole domain, otherwise you get wrong results.

I think you need to do both corrections: subtract the mean to f, and fix a single dof to 0.  
To compute the integral of f on the whole domain you need to use a partition of unity.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [January 28, 2026, 12:59pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/31 "2026-01-28T12:59:47Z")

</div>

Dear Prof. Bouchut,

Thank you for your suggestion. You are right — enforcing \int\_\Omega f = 0 is crucial. After subtracting the global mean of f, the solution becomes consistent.  
BTW, regarding `SetBC()`, which is correct, `tgv = -1` or `tgv = -2`?  
In particular, `tgv = -2` zero out both the row and the column associated with the constrained DoF, so that this DoF no longer contributes to the other DoFs?

Best regards,  
Haocheng Weng

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [January 28, 2026, 3:02pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/32 "2026-01-28T15:02:13Z")

</div>

tgv=-1 is always correct. Here the value -2 of tgv is also correct because the right-hand side (corresponding to the imposed dof) is zero. A good point for tgv=-2 is that it leads to a symmetric matrix, so that (depending on the method used for the linear system) it could give a more stable or faster resolution.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [January 29, 2026, 9:25am UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/33 "2026-01-29T09:25:24Z")

</div>

Dear Prof. Bouchut,

I noticed that once \int\_{\Omega} f = 0 is enforced, the solver seems able to return a solution even without fixing a DoF or explicitly specifying a null space. The solution is not unique, but it does satisfy the equation.

When only the pressure gradient is of interest, is it then unnecessary to use `SetBC()` or to explicitly fix the null space?

Best regards,  
Haocheng Weng

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [January 29, 2026, 3:12pm UTC](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351/34 "2026-01-29T15:12:13Z")

</div>

I think that it is necessary in principle to fix a dof in order to have an invertible matrix. However, when an iterative solver is used, it can converge to a solution even if there is no uniqueness, which is good for you. Hence you can skip fixing a dof as long as the solver does the job. If at some point it diverges, you know what you have to do!

[Previous page](https://community.freefem.org/t/solver-for-unsteady-n-s-equation/3351.md?page=1)