solving Poisson problem with pure Neumann conditions

[dahtah]

Is there any way of solving singular problems with amgcl? For example, a simple FD approximation of the Poisson problem with pure Neumann bc reads: [-1 1 .... ] [1 -2 1 .... ] [. 1 -2 1 ...] ... [ ..... 1 -1] and has a null space (functions with constant values). The standard CG algorithm converges on this problem when the rhs has zero mean. amgcl (with smoothed aggregation, spai0 relaxation) fails with:

Zero sum in skyline_lu factorization

which I'm guessing means the coarse system is singular. Do you know of a way of working around the problem?

[ddemidov]

which I'm guessing means the coarse system is singular. Do you know of a way of working around the problem?

That is correct.

You could change one of your equations to represent a fictitious Dirichlet condition (u(x)=0 at some arbitrary x). This should work in theory, because your solution is defined up to an additive constant, but in practice this could mean you get slower convergence rate. Also, most of your error would probably be located around the point you chose as an anchor. Unfortunately, I am not aware of any other options.

[ddemidov]

This tutorial suggests it should be possible to solve this kind of system by using smoothers and grid transfer operator with some postprocessing steps. I think this should be easy to add this functionality to amgcl. If you would you be interested in working on this, I could provide some pointers.

[dahtah]

I had a look at the tutorial and I'm not sure I understand the details. If we use CG all the way (incl. when solving coarse-grid problems), then it's enough to enforce the zero-mean constraint at each iteration. But if at some level of coarseness we start using direct solvers, we need to take into account the singularity of the equation. I imagine the LU decomposition could be replaced with a pseudo-inverse, or if a dense solver is used we can solve A+1*1^t instead of just A. I haven't used it, but as far as I understand the PETSC Krylov solver has the option of specifying a null space: http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetNullSpace.html I can't find any documentation on the algorithm, though.

[ddemidov]

You can supply a null-space vector for aggregation-type coarseners (see example here).

With the commit I pushed today you can also eliminate coarse direct solver entirely by setting max_levels and coarse_enough parameters to amgcl::amg so that number of equations on the coarsest level is always greater than coarse_enough. Then the coarsest level will apply a couple of smoother iterations instead of solving directly. Example:

with direct solver at coarse level

./examples/solver -n 64
coarsening:          smoothed_aggregation
relaxation:          spai0
Number of levels:    3
Operator complexity: 1.55
Grid complexity:     1.12

level     unknowns       nonzeros
---------------------------------
    0       262144        1810432 (64.45%)
    1        31868         955994 (34.04%)
    2          769          42421 ( 1.51%)

Iterations: 12
Error:      6.3746e-10

[Profile:          0.770 s] (100.00%)
[ self:            0.003 s] (  0.33%)
[  assembling:     0.022 s] (  2.79%)
[  setup:          0.250 s] ( 32.45%)
[  solve:          0.496 s] ( 64.42%)

with smoother at coarse level

./examples/solver -n 64 -p precond.coarse_enough=1 -p precond.max_levels=4
coarsening:          smoothed_aggregation
relaxation:          spai0
Number of levels:    4
Operator complexity: 1.55
Grid complexity:     1.12

level     unknowns       nonzeros
---------------------------------
    0       262144        1810432 (64.45%)
    1        31868         955994 (34.03%)
    2          769          42421 ( 1.51%)
    3           14            194 ( 0.01%)

Iterations: 13
Error:      6.41982e-09

[Profile:          0.809 s] (100.00%)
[ self:            0.003 s] (  0.36%)
[  assembling:     0.022 s] (  2.73%)
[  setup:          0.241 s] ( 29.76%)
[  solve:          0.543 s] ( 67.15%)

This may increase the number of iterations and the solution time. Unfortunately, I do not currently have a direct solver that could deal with a singular system.

EDIT: actually show output from master branch and not from a debug attempt.

[ddemidov]

You can also see if CG will work for you with a single-level smoother as preconditioner. The example of setting this up is here:

./examples/solver -n 64 -1
Using spai0 as preconditioner
  unknowns: 262144
  nonzeros: 1810432

Iterations: 100
Error:      2.66278e-08

[Profile:          0.998 s] (100.00%)
[ self:            0.003 s] (  0.29%)
[  assembling:     0.022 s] (  2.20%)
[  setup:          0.034 s] (  3.41%)
[  solve:          0.939 s] ( 94.11%)

[dahtah]

Great, thanks much! The first approach works, I still have to test the second one. Performance is pretty good on a toy Poisson problem.

[ddemidov]

Thanks for testing! b81db73 should allow to choose smoother over coarse direct solver more easily. It adds bool direct_coarse parameter. When that is set to false, a smoother is used instead of direct solver.