symmetric iterative linear solver malfunction?

The same works with UnitIntervalMesh(2) btw.

Works fine for me - check http://artax.karlin.mff.cuni.cz/~blecj6am/219418/
Altough solver reports different residual than assembler. But that's natural - solver calculates it somehow from matrix and rhs vector but and does not have information how much every dof contributes to integral by cell volume. Only assembler knows it.

Consiser also that there is only one dof not set by bcs at coords (0.5, 0.5).

Jan

This is the same output that I get. -- The final residual should be less than 1.0-10, not 0.125.

Add more unknowns to actually see the difference between the CG "solution" and the GMRES solution.

On Thu, 17 Jan 2013 13:50:59 -0000
Nico Schlömer <email address hidden> wrote:
> Question #219418 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/219418
>
> Nico Schlömer posted a new comment:
> This is the same output that I get. -- The final residual should be
> less than 1.0-10, not 0.125.
>
It can't be. In residual calculated by assembly of forms there is
included discretization error which is huge in this case (8 elements)!
But residual calculated by linear algebra backend is only
  norm(Ax-b)
where norm is some norm. In this simple case solver even achieved zero
residual after one iteration - it found exact solution of Ax=b -
altough Ax=b is approxamation (very rough!) to original problem.

What do you mean by "discretization error"? The between the exact solution of the full PDE problem and the discretized solution? That wouldn't play any role here.

I found out more, too: When explicitly assembling matrix and right-hand-side, and calling the Krylov solver on the linear problem, it all works out:

A = assemble(a, bcs=bcs)
b = assemble(L, bcs=bcs)
u = Function(V)
solve(A, u.vector(), b, 'cg')

This definitely seems like a bug in the solve(a == L, ...).

> What do you mean by "discretization error"? The between the exact
> solution of the full PDE problem and the discretized solution?
Yes.

> That
> wouldn't play any role here.
I don't agree. Doing
> # Compute the residual norm.
> r = zero() * dx(0)
> for i in [0,1,2]:
> r += inner(K[i]*grad(u), grad(v)) * dx(i) \
> - f*v*dx(i)
> R = assemble(r, bcs=bcs)
you calculate residual in original continuous problem with u discrete
solution. That's ensured because ffc uses qaudrature order ensuring
exact qudrature in this case.

Again this solution http://artax.karlin.mff.cuni.cz/~blecj6am/219418/
is exact solution to discrete problem (zero linear algebra residual) but
with a huge discretization error (non-zero residual by assembling
forms)! Very simple...

> This definitely seems like a bug in the solve(a == L, ...).
>
It does not for me.

>> r += inner(K[i]*grad(u), grad(v)) * dx(i) \
>> - f*v*dx(i)
> you calculate residual in original continuous problem with u discrete solution.

Hm, I wouldn't think so. u and v are already in the discretized space V, and assemble() iterates over them just like the assemble() step in the creation of the linear system does.
Besides, this residual *is* of the order of 1.0e-18 when using GMRES or when K is constant throughout the domain.

Also, the solution that you show is *not* the solution of the discretized system. Add more points to the discretization and it will become clearer.

> >> r += inner(K[i]*grad(u), grad(v)) * dx(i) \
> >> - f*v*dx(i)
> > you calculate residual in original continuous problem with u
> > discrete solution.
>
> Hm, I wouldn't think so. u and v are already in the discretized space
> V, and assemble() iterates over them just like the assemble() step in
> the creation of the linear system does. Besides, this residual *is*
> of the order of 1.0e-18 when using GMRES or when K is constant
> throughout the domain.
In fact you're right. Sorry, my mistake.

But I think I've revealed a problem. First of all you've got wrong
definition of Left and Right from a numerical point of view. Try
    x[0] <= 0.5 + DOLFIN_EPS
    0.5 - DOLFIN_EPS <= x[0]

But importantly it seems being a bug. It seems that if you set
'symmetric': True, solve method calls assemble_system() instead of
assemble(). And if you try

    A, b = assemble_system(a, L, bcs=bcs)
    solve(A, u.vector(), b, 'cg', 'none')

you get same wrong beviour as with solve(a==L, ...{'symmetric': True}).
But with

    A, b = assemble_system(a, L, bcs=bcs, cell_domains=domains)
    solve(A, u.vector(), b, 'cg', 'none')

it works. So it seems that assemble_system() method ignores
cell_domains stored in Form and only uses explicitly specified
cell_domains. If you confirm it we should report a bug.

Interesting is that cg solver works also on asymmetric matrix as you
suggested. Try

    A = assemble(a, bcs=bcs)
    b = assemble(L, bcs=bcs)
    solve(A, u.vector(), b, 'cg', 'none')

assemble() method does not ignore cell_domains attached to a and L.

Jan

Thanks for looking into this!

> If you confirm it we should report a bug.

I did just that.

> Interesting is that cg solver works also on asymmetric matrix as you suggested.

This is only for special cases. The case here is:

The matrix has (up to reordering) the form

( I 0 )
( A B )

where B is spd *and* the initial guess (zeros()) fulfills the boundary condition. As a consequence, the initial residual of the boundary part of u is 0 and will never play a role throughout the iteration. Morale: Adapt your initial guesses! :)