Dolfin GMRES vs. CBC.Block LGMRES

I do not have any ideas off-hand. I have used LGMRES without seing this
problem before. Which dolfin and cbc.block version are you running? Also,
if you could provide a runnable example that shows this error (complete,
not necessarily minimal) then I can have a look. Send to <email address hidden>.

On 28 January 2013 12:06, Claas Abert
<email address hidden>wrote:

> New question #220333 on CBC.Block:
> https://answers.launchpad.net/cbc.block/+question/220333
>
> This question is maybe way too general, but perhaps you can give me some
> hints:
>
> I want to solve a nonsymmetric problem with CBC.Block.
> The standard dolfin solve method with Trilinos backend solves the system
> in 5 iterations:
>
> A, b = assemble_system(a, L)
> solve(A, x.vector(), b, "gmres", "ilu")
>
> However trying to solve the system with CBC.Block fails:
> Ap = ILU(A)
> Ainv = LGMRES(A, precond=Ap)
> x = Ainv * b
>
> with the following:
>
> File "/usr/local/lib/python2.7/dist-packages/block/iterative/lgmres.py",
> line 190, in lgmres
> y, resids, rank, s = lstsq(hess, e1)
> File "/usr/lib/python2.7/dist-packages/scipy/linalg/basic.py", line 409,
> in lstsq
> a1, b1 = map(asarray_chkfinite, (a, b))
> File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line
> 590, in asarray_chkfinite
> "array must not contain infs or NaNs")
> ValueError: array must not contain infs or NaNs
>
> I also tried BiCGStab as a solver (3 iterations with builtin Trilinos
> solver, NaN residuum in case of the CBC.Block implementation).
>
> As far as I understand, both the Trilinos GMRES and the CBC.Block LGMRES
> use the IFPACK ILU preconditioner. So does this mean the only difference
> between the two libraries is the implementation of the iterative methods?
>
> I want to use the CBC.Block solver, because I want to implement a Schur
> complement method where I have to solve another linear system per GMRES
> step.
>
> Do you have any hints what I could do to make things work with CBC.Block?
> I know this question is very general. I will gladly provide more
> information if needed.
>
> Thanks in advance, Claas
>
> --
> You received this question notification because you are a member of
> cbc.block maintainers, which is an answer contact for CBC.Block.
>
> --
> Mailing list: https://launchpad.net/~cbc.block
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~cbc.block
> More help : https://help.launchpad.net/ListHelp
>

I can also look at it, so please post to cbc.block if the code
is not unreasonably large

Kent

On 28 January 2013 12:15, Joachim Haga <<email address hidden>
> wrote:

> Question #220333 on CBC.Block changed:
> https://answers.launchpad.net/cbc.block/+question/220333
>
> Status: Open => Answered
>
> Joachim Haga proposed the following answer:
> I do not have any ideas off-hand. I have used LGMRES without seing this
> problem before. Which dolfin and cbc.block version are you running? Also,
> if you could provide a runnable example that shows this error (complete,
> not necessarily minimal) then I can have a look. Send to <email address hidden>.
>
>
> On 28 January 2013 12:06, Claas Abert
> <email address hidden>wrote:
>
> > New question #220333 on CBC.Block:
> > https://answers.launchpad.net/cbc.block/+question/220333
> >
> > This question is maybe way too general, but perhaps you can give me some
> > hints:
> >
> > I want to solve a nonsymmetric problem with CBC.Block.
> > The standard dolfin solve method with Trilinos backend solves the system
> > in 5 iterations:
> >
> > A, b = assemble_system(a, L)
> > solve(A, x.vector(), b, "gmres", "ilu")
> >
> > However trying to solve the system with CBC.Block fails:
> > Ap = ILU(A)
> > Ainv = LGMRES(A, precond=Ap)
> > x = Ainv * b
> >
> > with the following:
> >
> > File "/usr/local/lib/python2.7/dist-packages/block/iterative/lgmres.py",
> > line 190, in lgmres
> > y, resids, rank, s = lstsq(hess, e1)
> > File "/usr/lib/python2.7/dist-packages/scipy/linalg/basic.py", line 409,
> > in lstsq
> > a1, b1 = map(asarray_chkfinite, (a, b))
> > File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line
> > 590, in asarray_chkfinite
> > "array must not contain infs or NaNs")
> > ValueError: array must not contain infs or NaNs
> >
> > I also tried BiCGStab as a solver (3 iterations with builtin Trilinos
> > solver, NaN residuum in case of the CBC.Block implementation).
> >
> > As far as I understand, both the Trilinos GMRES and the CBC.Block LGMRES
> > use the IFPACK ILU preconditioner. So does this mean the only difference
> > between the two libraries is the implementation of the iterative methods?
> >
> > I want to use the CBC.Block solver, because I want to implement a Schur
> > complement method where I have to solve another linear system per GMRES
> > step.
> >
> > Do you have any hints what I could do to make things work with CBC.Block?
> > I know this question is very general. I will gladly provide more
> > information if needed.
> >
> > Thanks in advance, Claas
> >
> > --
> > You received this question notification because you are a member of
> > cbc.block maintainers, which is an answer contact for CBC.Block.
> >
> > --
> > Mailing list: https://launchpad.net/~cbc.block
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~cbc.block
> > More help : https://help.launchpad.net/ListHelp
> >
>
> --
> You received this question notification because you are a member of
> cbc.block maintainers, which is an answer contact for CBC.Block.
>
> --
> Mailing list: https://launchpad.net/~cbc.block
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~cbc.block
> More help : https://help.launchpad.net/ListHelp
>

Thanks Joachim Haga and Kent-Andre Mardal for the quick response.
I'm using the dolfin-dev rev 7250 (from the beginning of 2013) and the most recent dev-version of cbc.block (rev 121).
Here is a (not very minimal) example showing the problem.

The problem itself is a saddle point problem and I know that I should probably deal with this block structure properly using CBC.Block. However this is just the starting point for a bigger block system and I wanted to keep things simple in the beginning...

from dolfin import *
from block import *
from block.iterative import *
from block.algebraic.trilinos import *

# Set backend to Epetra
parameters['linear_algebra_backend'] = 'Epetra'

# Mesh and Function Spaces
mesh = BoxMesh(-50, -50, -5, 50, 50, 5, 10, 10, 1)
VV = VectorFunctionSpace(mesh, "CG", 1)
VS = FunctionSpace(mesh, "CG", 1)
V = VV * VS

# Test and Trial Functions
(v, sigma) = TrialFunctions(V)
(w2, w3) = TestFunctions(V)

# Define some constants and coeffecients
arg = "sqrt((3.141592*(x[0]/1e1))*(3.141592*(x[0]/1e1)))"
m_expr = Expression(("cos(%s)" % arg, "sin(%s)" % arg, "0.0"))
m = interpolate(m_expr, VV)

f_ex = -5.71823816178e+12
alpha = 0.1
dt = 1e-12

# Define weak forms
a = alpha * dot(v, w2) * dx + dot(cross(m, v), w2) * dx
a += - 0.5 * dt * f_ex * Dx(v[i],j) * Dx(w2[i],j) * dx
a += sigma * inner(m, w2) * dx
a += inner(m, v) * w3 * dx

L = f_ex * Dx(m[i],j) * Dx(w2[i],j) * dx

# Assemble
A = assemble(a)
b = assemble(L)

# Solve with Trilinos solvers
x = Function(V)
solve(A, x.vector(), b, "gmres", "ilu")
solve(A, x.vector(), b, "bicgstab", "ilu")

# Solve with CBC.Block solvers
Ap = ILU(A)
Ainv = LGMRES(A, precond=Ap)
#Ainv = BiCGStab(A, precond=Ap)
x = Ainv * b

I have looked briefly. The convergence of builtin Trilinos solvers is illusory; try

# Solve with Trilinos solvers
x = Function(V)
solve(A, x.vector(), b, "gmres", "ilu")
plot(x[0])
solve(A, x.vector(), b, "bicgstab", "ilu")
plot(x[0])
solve(A, x.vector(), b, "bicgstab", "ilu")
plot(x[0])
interactive()

--- you'll see three very different "solutions". I suspect that the problem is that the full matrix A is indefinite; the preconditioner should probably be positive. The built-in solver probably bails out when encountering division-by-zero (which generates NaN/Inf), while cbc.block doesn't check for this and hence throws an exception.

Oops, that last solve() should be just

solve(A, x.vector(), b)

to use a direct solver.

Anyway, it looks like simple doesn't work here. Go blockwise ;)

Thanks a lot Joachim!
I'm sorry that I didn't check that before posting.
I will go for the blocks now :)