ddemidov / amgcl

C++ library for solving large sparse linear systems with algebraic multigrid method
http://amgcl.readthedocs.org/
MIT License
728 stars 111 forks source link

Convergence issue when directly using AMG as a linear solver #97

Closed jiongchen closed 3 years ago

jiongchen commented 5 years ago

Dear Denis,

I really would like to thank you for open sourcing such a well designed and super flexible framework to play with AMG:)

I was working on a project on multigrid recently, and tried to make some comparisons with the standard AMG algorithms. Instead of using AMG as a preconditioner for an iterative linear solver, I tried to solve the Poisson equation with amg class directly (by calling apply method). But somehow I found that it is hard to decrease the error but actually increase it a bit instead, no matter ruge stuben or smoothed aggregation coarsening strategies is applied. Some parameters are set as (using damped jacobi as smoother):

npre = npost = 3;
relax.damping = 0.3;
max_levels = 3;
pre_cycles = 100;

Since I am not an expert on AMG, I cannot figure out the reason behind. But to my understanding, the homogeneous Laplacian matrix should not be challenging for AMG. So is it a normal phenomena? Or just caused by the inappropriate setting of the parameters?

ddemidov commented 5 years ago

smg::apply() is meant to be used as preconditioning step. If you look at the code, it starts by zeroing out the x vector, and then calls amg::cycle.

For your use case you should call amg::cycle directly:

https://github.com/ddemidov/amgcl/blob/91348b025144b61a2d3b7417988373d0dc8c8d00/amgcl/amg.hpp#L255-L263

jiongchen commented 5 years ago

For the Poisson equation given in the example, I initialize the solution vector by zero, and iteratively call amg::cycle for about 100 times. Still, the residual norm will increase instead of decrease. Have you ever encountered such a problem?

ddemidov commented 5 years ago

Here is the solution that works for me: https://gist.github.com/b05ce89393718647c2b1d57fa2ad0043

The output I get from this is:

Number of levels:    3
Operator complexity: 1.59
Grid complexity:     1.14
Memory footprint:    11.54 M

level     unknowns       nonzeros      memory
---------------------------------------------
    0        32768         223232      8.48 M (62.80%)
    1         4192         114356      2.49 M (32.17%)
    2          346          17898    579.04 K ( 5.03%)

Iterations: 20
Error     : 6.29241e-07

[Profile:     1.006 s] (100.00%)
[ self:       0.001 s] (  0.13%)
[  setup:     0.133 s] ( 13.22%)
[  solve:     0.872 s] ( 86.65%)

When solving the same problem with BiCGStab preconditioned by AMG, I get:

Solver
======
Type:             BiCGStab
Unknowns:         32768
Memory footprint: 1.75 M

Preconditioner
==============
Number of levels:    3
Operator complexity: 1.59
Grid complexity:     1.14
Memory footprint:    11.54 M

level     unknowns       nonzeros      memory
---------------------------------------------
    0        32768         223232      8.48 M (62.80%)
    1         4192         114356      2.49 M (32.17%)
    2          346          17898    579.04 K ( 5.03%)

Iterations: 6
Error:      2.73764e-07

[Profile:          0.765 s] (100.00%)
[ self:            0.005 s] (  0.64%)
[  assembling:     0.002 s] (  0.29%)
[  setup:          0.148 s] ( 19.39%)
[  solve:          0.609 s] ( 79.68%)

(6 iterations instead of 20).

I am using default parameters here. In your initial comment you set damping factor to 0.3. That could be too low (I think something like 0.72 should work better).

jiongchen commented 5 years ago

I will try to run the given code. Maybe I made some stupid errors in my own example. It's very kind of you and I really appreciate for your help! Thank you so much!

ddemidov commented 3 years ago

There is a new Richardson iterative solver, which should be equivalent to applying the preconditioner in a loop until convergence.