[petsc-users] Nullspace issue

[Thomas Witkowski]

I've some trouble to setup the nullspace for my equation correctly. It
is an implicit coupled Navier-Stokes / Cahn-Hilliard equation (for two
phase flow simulation) in 2D. To simulate an "infinite" domain, the
upper and lower boundaries are periodic for all 5 components. The
pressure and the two variables from the Cahn-Hilliard equation are
null-flux (Neumann) on the other two boundaries, the velocity is set to
0 (Dirichlet) on them. I would expect to have a non-empty nullspace.
When using richardson/lu, I see the following behavior:

1 KSP preconditioned resid norm 3.130958085604e+04 true resid norm
1.339173230263e-11 ||r(i)||/||b|| 2.249831332702e-16
2 KSP preconditioned resid norm 5.987054624984e+04 true resid norm
8.900173028027e-12 ||r(i)||/||b|| 1.495242564025e-16
3 KSP preconditioned resid norm 8.159870664335e+04 true resid norm
7.429941832854e-12 ||r(i)||/||b|| 1.248241493927e-16
4 KSP preconditioned resid norm 7.087492912608e+04 true resid norm
7.066574102540e-12 ||r(i)||/||b|| 1.187195164260e-16
5 KSP preconditioned resid norm 1.075156685382e+05 true resid norm
6.963290776604e-12 ||r(i)||/||b|| 1.169843408895e-16
6 KSP preconditioned resid norm 9.260088828115e+04 true resid norm
6.996993371924e-12 ||r(i)||/||b|| 1.175505496012e-16

So, in principle the system is solvable, but the preconditioned norm is
wrong. I created a nullspace basis that is 1 on all pressure unknowns
and 0 on all of the other four components. I can easily verify that this
vector is at least an element of the nullspace, as the norm of the
vector resulting from MatMult with the assembled matrix and the
nullspace basis vector is in the order of 1e-16. I set the nullspace to
the system matrix, and rerun with richardson/lu:

0 KSP preconditioned resid norm 6.025198436485e+09 true resid norm
5.952327229149e+04 ||r(i)||/||b|| 1.000000000000e+00
1 KSP preconditioned resid norm 3.436455497312e+09 true resid norm
9.795164213312e-09 ||r(i)||/||b|| 1.645602440226e-13
2 KSP preconditioned resid norm 3.365625775916e+09 true resid norm
7.488946142386e-09 ||r(i)||/||b|| 1.258154307396e-13
3 KSP preconditioned resid norm 2.543011385543e+09 true resid norm
5.719959416026e-09 ||r(i)||/||b|| 9.609618550566e-14
4 KSP preconditioned resid norm 4.406639844892e+08 true resid norm
6.417414971975e-09 ||r(i)||/||b|| 1.078135446007e-13
5 KSP preconditioned resid norm 2.741735667402e+09 true resid norm
5.646995693759e-09 ||r(i)||/||b|| 9.487038390809e-14
6 KSP preconditioned resid norm 2.427193596206e+09 true resid norm
1.327430895972e-08 ||r(i)||/||b|| 2.230104033043e-13
7 KSP preconditioned resid norm 1.603502416252e+09 true resid norm
1.338139212677e-08 ||r(i)||/||b|| 2.248094167478e-13
8 KSP preconditioned resid norm 1.646847388237e+09 true resid norm
1.287770565938e-08 ||r(i)||/||b|| 2.163474077218e-13
9 KSP preconditioned resid norm 2.136833312630e+09 true resid norm
1.261183165000e-08 ||r(i)||/||b|| 2.118806840498e-13
10 KSP preconditioned resid norm 1.405267823532e+09 true resid norm
2.499790590170e-08 ||r(i)||/||b|| 4.199686095765e-13

The system at all is not solvable anymore. What could be the problem?
Some technical issue I forgot about? Or my be the nullspace be wrong?
Thanks for any hint!

Thomas

[Jed Brown]

A known null space doesn't really "fix" LU for a singular problem [1]. You
need a preconditioner that is *stable*. Also, it sounds like you may have
additional vectors in the null space coming from the Cahn-Hilliard
variables.

[1] One could create a new matrix with space for Lagrange multipliers
enforcing the null space and perform the factorization there. But doing
this automatically isn't implemented in the library so you'd have to write
it.

On Tue, Jun 12, 2012 at 1:32 AM, Thomas Witkowski <

I've some trouble to setup the nullspace for my equation correctly. It is
an implicit coupled Navier-Stokes / Cahn-Hilliard equation (for two phase
flow simulation) in 2D. To simulate an "infinite" domain, the upper and
lower boundaries are periodic for all 5 components. The pressure and the
two variables from the Cahn-Hilliard equation are null-flux (Neumann) on
the other two boundaries, the velocity is set to 0 (Dirichlet) on them. I
would expect to have a non-empty nullspace. When using richardson/lu, I see
1 KSP preconditioned resid norm 3.130958085604e+04 true resid norm
1.339173230263e-11 ||r(i)||/||b|| 2.249831332702e-16
2 KSP preconditioned resid norm 5.987054624984e+04 true resid norm
8.900173028027e-12 ||r(i)||/||b|| 1.495242564025e-16
3 KSP preconditioned resid norm 8.159870664335e+04 true resid norm
7.429941832854e-12 ||r(i)||/||b|| 1.248241493927e-16
4 KSP preconditioned resid norm 7.087492912608e+04 true resid norm
7.066574102540e-12 ||r(i)||/||b|| 1.187195164260e-16
5 KSP preconditioned resid norm 1.075156685382e+05 true resid norm
6.963290776604e-12 ||r(i)||/||b|| 1.169843408895e-16
6 KSP preconditioned resid norm 9.260088828115e+04 true resid norm
6.996993371924e-12 ||r(i)||/||b|| 1.175505496012e-16
So, in principle the system is solvable, but the preconditioned norm is
wrong. I created a nullspace basis that is 1 on all pressure unknowns and 0
on all of the other four components. I can easily verify that this vector
is at least an element of the nullspace, as the norm of the vector
resulting from MatMult with the assembled matrix and the nullspace basis
vector is in the order of 1e-16. I set the nullspace to the system matrix,
0 KSP preconditioned resid norm 6.025198436485e+09 true resid norm
5.952327229149e+04 ||r(i)||/||b|| 1.000000000000e+00
1 KSP preconditioned resid norm 3.436455497312e+09 true resid norm
9.795164213312e-09 ||r(i)||/||b|| 1.645602440226e-13
2 KSP preconditioned resid norm 3.365625775916e+09 true resid norm
7.488946142386e-09 ||r(i)||/||b|| 1.258154307396e-13
3 KSP preconditioned resid norm 2.543011385543e+09 true resid norm
5.719959416026e-09 ||r(i)||/||b|| 9.609618550566e-14
4 KSP preconditioned resid norm 4.406639844892e+08 true resid norm
6.417414971975e-09 ||r(i)||/||b|| 1.078135446007e-13
5 KSP preconditioned resid norm 2.741735667402e+09 true resid norm
5.646995693759e-09 ||r(i)||/||b|| 9.487038390809e-14
6 KSP preconditioned resid norm 2.427193596206e+09 true resid norm
1.327430895972e-08 ||r(i)||/||b|| 2.230104033043e-13
7 KSP preconditioned resid norm 1.603502416252e+09 true resid norm
1.338139212677e-08 ||r(i)||/||b|| 2.248094167478e-13
8 KSP preconditioned resid norm 1.646847388237e+09 true resid norm
1.287770565938e-08 ||r(i)||/||b|| 2.163474077218e-13
9 KSP preconditioned resid norm 2.136833312630e+09 true resid norm
1.261183165000e-08 ||r(i)||/||b|| 2.118806840498e-13
10 KSP preconditioned resid norm 1.405267823532e+09 true resid norm
2.499790590170e-08 ||r(i)||/||b|| 4.199686095765e-13
The system at all is not solvable anymore. What could be the problem? Some
technical issue I forgot about? Or my be the nullspace be wrong? Thanks for
any hint!
Thomas

[Thomas Witkowski]

There should be no null space from the Cahn-Hilliard equation. Is there
some black-box preconditioner that does not relay on LU factorization at
some point? I know that black-box approaches are mostly not efficient,
but I would have something I can work with.

Thomas

Post by Jed Brown
A known null space doesn't really "fix" LU for a singular problem [1].
You need a preconditioner that is *stable*. Also, it sounds like you
may have additional vectors in the null space coming from the
Cahn-Hilliard variables.
[1] One could create a new matrix with space for Lagrange multipliers
enforcing the null space and perform the factorization there. But
doing this automatically isn't implemented in the library so you'd
have to write it.
On Tue, Jun 12, 2012 at 1:32 AM, Thomas Witkowski
I've some trouble to setup the nullspace for my equation
correctly. It is an implicit coupled Navier-Stokes / Cahn-Hilliard
equation (for two phase flow simulation) in 2D. To simulate an
"infinite" domain, the upper and lower boundaries are periodic for
all 5 components. The pressure and the two variables from the
Cahn-Hilliard equation are null-flux (Neumann) on the other two
boundaries, the velocity is set to 0 (Dirichlet) on them. I would
expect to have a non-empty nullspace. When using richardson/lu, I
1 KSP preconditioned resid norm 3.130958085604e+04 true resid
norm 1.339173230263e-11 ||r(i)||/||b|| 2.249831332702e-16
2 KSP preconditioned resid norm 5.987054624984e+04 true resid
norm 8.900173028027e-12 ||r(i)||/||b|| 1.495242564025e-16
3 KSP preconditioned resid norm 8.159870664335e+04 true resid
norm 7.429941832854e-12 ||r(i)||/||b|| 1.248241493927e-16
4 KSP preconditioned resid norm 7.087492912608e+04 true resid
norm 7.066574102540e-12 ||r(i)||/||b|| 1.187195164260e-16
5 KSP preconditioned resid norm 1.075156685382e+05 true resid
norm 6.963290776604e-12 ||r(i)||/||b|| 1.169843408895e-16
6 KSP preconditioned resid norm 9.260088828115e+04 true resid
norm 6.996993371924e-12 ||r(i)||/||b|| 1.175505496012e-16
So, in principle the system is solvable, but the preconditioned
norm is wrong. I created a nullspace basis that is 1 on all
pressure unknowns and 0 on all of the other four components. I can
easily verify that this vector is at least an element of the
nullspace, as the norm of the vector resulting from MatMult with
the assembled matrix and the nullspace basis vector is in the
order of 1e-16. I set the nullspace to the system matrix, and
0 KSP preconditioned resid norm 6.025198436485e+09 true resid
norm 5.952327229149e+04 ||r(i)||/||b|| 1.000000000000e+00
1 KSP preconditioned resid norm 3.436455497312e+09 true resid
norm 9.795164213312e-09 ||r(i)||/||b|| 1.645602440226e-13
2 KSP preconditioned resid norm 3.365625775916e+09 true resid
norm 7.488946142386e-09 ||r(i)||/||b|| 1.258154307396e-13
3 KSP preconditioned resid norm 2.543011385543e+09 true resid
norm 5.719959416026e-09 ||r(i)||/||b|| 9.609618550566e-14
4 KSP preconditioned resid norm 4.406639844892e+08 true resid
norm 6.417414971975e-09 ||r(i)||/||b|| 1.078135446007e-13
5 KSP preconditioned resid norm 2.741735667402e+09 true resid
norm 5.646995693759e-09 ||r(i)||/||b|| 9.487038390809e-14
6 KSP preconditioned resid norm 2.427193596206e+09 true resid
norm 1.327430895972e-08 ||r(i)||/||b|| 2.230104033043e-13
7 KSP preconditioned resid norm 1.603502416252e+09 true resid
norm 1.338139212677e-08 ||r(i)||/||b|| 2.248094167478e-13
8 KSP preconditioned resid norm 1.646847388237e+09 true resid
norm 1.287770565938e-08 ||r(i)||/||b|| 2.163474077218e-13
9 KSP preconditioned resid norm 2.136833312630e+09 true resid
norm 1.261183165000e-08 ||r(i)||/||b|| 2.118806840498e-13
10 KSP preconditioned resid norm 1.405267823532e+09 true resid
norm 2.499790590170e-08 ||r(i)||/||b|| 4.199686095765e-13
The system at all is not solvable anymore. What could be the
problem? Some technical issue I forgot about? Or my be the
nullspace be wrong? Thanks for any hint!
Thomas

[Jed Brown]

On Tue, Jun 12, 2012 at 7:56 AM, Thomas Witkowski <

Post by Thomas Witkowski
There should be no null space from the Cahn-Hilliard equation.

You said all those boundary conditions are either Neumann or periodic. I
guess it couples to the fluid variables without any null space?

Post by Thomas Witkowski
Is there some black-box preconditioner that does not relay on LU
factorization at some point? I know that black-box approaches are mostly
not efficient, but I would have something I can work with.

The SVD always works and will tell you about a null space, but of course
it's very expensive.

[Thomas Witkowski]

Post by Jed Brown
On Tue, Jun 12, 2012 at 7:56 AM, Thomas Witkowski
There should be no null space from the Cahn-Hilliard equation.
You said all those boundary conditions are either Neumann or periodic.
I guess it couples to the fluid variables without any null space?

yes.

Post by Jed Brown
Is there some black-box preconditioner that does not relay on LU
factorization at some point? I know that black-box approaches are
mostly not efficient, but I would have something I can work with.
The SVD always works and will tell you about a null space, but of
course it's very expensive.

So assume I have a basis for the null space of the system that should be
solved. Is there any block-box solver/preconditioner approach that does
not make use of (I)LU factorization at any point?

Thomas

[Jed Brown]

SVD works. Black box solvers aren't really.

You can try sparse approximate inverse, but good luck. ;-D

If you specify a "near null space", you can try smoothed aggregation.

Or just run domain decomposition, the decomposition generally breaks to
null space so that factorization can be used on subdomains.

Post by Jed Brown
On Tue, Jun 12, 2012 at 7:56 AM, Thomas Witkowski <

Post by Thomas Witkowski
There should be no null space from the Cahn-Hilliard equation.

You said all those boundary conditions are either Neumann or periodic. I
guess it couples to the fluid variables without any null space?
yes.

Post by Thomas Witkowski
Is there some black-box preconditioner that does not relay on LU
factorization at some point? I know that black-box approaches are mostly
not efficient, but I would have something I can work with.

The SVD always works and will tell you about a null space, but of course
it's very expensive.
So assume I have a basis for the null space of the system that should be
solved. Is there any block-box solver/preconditioner approach that does not
make use of (I)LU factorization at any point?
Thomas

[Barry Smith]

Thomas,

Since you are solving a coupled system of equations you should start with PCFIELDSPLIT. This allows you to build a preconditioner by combining solvers for separate components of the PDE. It can even be applied recursively to first separate the Navier-Stokes equations from the Cahn-Hillard and then to separate the parts of the Navier-Stokes.

It may take a little thought for how you want the separation done and what solvers to use on what parts, but it is the only way to get an efficient solver for such a coupled system.

Barry

Jed, WTF are you talking about SVDs and stuff?

Post by Thomas Witkowski
There should be no null space from the Cahn-Hilliard equation.
You said all those boundary conditions are either Neumann or periodic. I guess it couples to the fluid variables without any null space?

yes.

Post by Thomas Witkowski
Is there some black-box preconditioner that does not relay on LU factorization at some point? I know that black-box approaches are mostly not efficient, but I would have something I can work with.
The SVD always works and will tell you about a null space, but of course it's very expensive.

So assume I have a basis for the null space of the system that should be solved. Is there any block-box solver/preconditioner approach that does not make use of (I)LU factorization at any point?
Thomas

[Jed Brown]

Post by Barry Smith
Jed, WTF are you talking about SVDs and stuff?

It's "black box" and works, if you are sufficiently unambitious and have
endless time to wait.

[Thomas Witkowski]

[Barry Smith]

From the three available block algorithms of the fieldsplit preconditioner, I could use only block jacobi and gauss seidel ones. For the Schur type approach I have no idea what to do with the Schur complement in the case block one is the Navier Stokes equation and block two is the Cahn Hilliard equation. Is there any assumption (that could also be verified in practice) on the system that these block preconditioners will work? What is your experience? If there are good preconditioners for the two blocks available, how good will pcfieldsplit preconditioner work?

G-S and Shur complement type fieldsplit PCs almost always work right so long as you have decent preconditioners for the original block(s). Often you don't need much of anything to precondition the Schur complement.

As I said before you can construct a preconditioner for the N.S. block by again using a PCFIELDSPLIT on that (recursive use of PCs).

Note that once you have defined your blocks (in the simple case if you have all your degrees of freedom collocated (no staggered grid) just set the Vec and Mat block size; otherwise you need to construct an IS that defines each block) then experimenting is easy and can be done at the command line without changing the code.

Barry

Thomas

Post by Barry Smith
Thomas,
Since you are solving a coupled system of equations you should start with PCFIELDSPLIT. This allows you to build a preconditioner by combining solvers for separate components of the PDE. It can even be applied recursively to first separate the Navier-Stokes equations from the Cahn-Hillard and then to separate the parts of the Navier-Stokes.
It may take a little thought for how you want the separation done and what solvers to use on what parts, but it is the only way to get an efficient solver for such a coupled system.
Barry
Jed, WTF are you talking about SVDs and stuff?

Post by Thomas Witkowski
There should be no null space from the Cahn-Hilliard equation.
You said all those boundary conditions are either Neumann or periodic. I guess it couples to the fluid variables without any null space?

yes.

Post by Thomas Witkowski
Is there some black-box preconditioner that does not relay on LU factorization at some point? I know that black-box approaches are mostly not efficient, but I would have something I can work with.
The SVD always works and will tell you about a null space, but of course it's very expensive.

So assume I have a basis for the null space of the system that should be solved. Is there any block-box solver/preconditioner approach that does not make use of (I)LU factorization at any point?
Thomas

[Jed Brown]

Post by Barry Smith
G-S and Shur complement type fieldsplit PCs almost always work right so
long as you have decent preconditioners for the original block(s). Often
you don't need much of anything to precondition the Schur complement.

Thomas, note that you'll need Schur or some surrogate preconditioner to
deal with incompressibility. It'll be worth understanding fieldsplit
variants while working with Navier-Stokes and Cahn-Hilliard separately.
Once you understand how to "drive" fieldsplit on those systems separately,
you can put them together.

Post by Barry Smith
As I said before you can construct a preconditioner for the N.S. block
by again using a PCFIELDSPLIT on that (recursive use of PCs).
Note that once you have defined your blocks (in the simple case if you
have all your degrees of freedom collocated (no staggered grid) just set
the Vec and Mat block size; otherwise you need to construct an IS that
defines each block) then experimenting is easy and can be done at the
command line without changing the code.

[Barry Smith]

Post by Barry Smith
G-S and Shur complement type fieldsplit PCs almost always work right so long as you have decent preconditioners for the original block(s). Often you don't need much of anything to precondition the Schur complement.
Thomas, note that you'll need Schur or some surrogate preconditioner to deal with incompressibility. It'll be worth understanding fieldsplit variants while working with Navier-Stokes and Cahn-Hilliard separately. Once you understand how to "drive" fieldsplit on those systems separately, you can put them together.

Very good point Jed. Until you can solve N.S. efficiently (for example using PCFIELDSPLIT etc) you really cannot put them together. But that exercise is not "wasted" effort because you will understand fieldsplit and Schur preconditioners and then extending to C.H will be much easier. PCFIELDSPLIT is very much a divide and conquer algorithmic approach, for for divide and conquer to work you need to know how to conquer the pieces.

Barry

Post by Barry Smith
As I said before you can construct a preconditioner for the N.S. block by again using a PCFIELDSPLIT on that (recursive use of PCs).
Note that once you have defined your blocks (in the simple case if you have all your degrees of freedom collocated (no staggered grid) just set the Vec and Mat block size; otherwise you need to construct an IS that defines each block) then experimenting is easy and can be done at the command line without changing the code.

[Thomas Witkowski]

Barry and Jed: Thank you very much for all the hints. Now I have a clue
where I have to start ...

Thomas

Post by Barry Smith

Post by Barry Smith
G-S and Shur complement type fieldsplit PCs almost always work right so long as you have decent preconditioners for the original block(s). Often you don't need much of anything to precondition the Schur complement.
Thomas, note that you'll need Schur or some surrogate preconditioner to deal with incompressibility. It'll be worth understanding fieldsplit variants while working with Navier-Stokes and Cahn-Hilliard separately. Once you understand how to "drive" fieldsplit on those systems separately, you can put them together.

Very good point Jed. Until you can solve N.S. efficiently (for example using PCFIELDSPLIT etc) you really cannot put them together. But that exercise is not "wasted" effort because you will understand fieldsplit and Schur preconditioners and then extending to C.H will be much easier. PCFIELDSPLIT is very much a divide and conquer algorithmic approach, for for divide and conquer to work you need to know how to conquer the pieces.
Barry

Post by Barry Smith
As I said before you can construct a preconditioner for the N.S. block by again using a PCFIELDSPLIT on that (recursive use of PCs).
Note that once you have defined your blocks (in the simple case if you have all your degrees of freedom collocated (no staggered grid) just set the Vec and Mat block size; otherwise you need to construct an IS that defines each block) then experimenting is easy and can be done at the command line without changing the code.