Here's a complete example:
=============== *snip* ===============
from dolfin import *
mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 1)
f = Expression('sin(x[0]) * sin(x[1])')
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
L = f * v * dx
A = assemble(a)
b = assemble(L)
# Make sure the system is consistent, i.e., the right-hand side is in the image
# of A. This is equivalent to the rhs not having a component in the kernel of
# A (since A is symmetric).
e = Function(V)
e.interpolate(Constant(1.0))
evec = e.vector()
evec /= norm(evec)
alpha = b.inner(evec)
b -= alpha * evec
x = Function(V)
prec = PETScPreconditioner('hypre_amg')
prec.parameters['hypre']['BoomerAMG']['relax_type_coarse'] = 'jacobi'
solver = PETScKrylovSolver('cg')
solver.parameters['absolute_tolerance'] = 0.0
solver.parameters['relative_tolerance'] = 1.0e-10
solver.parameters['maximum_iterations'] = 100
solver.parameters['monitor_convergence'] = True
# Create solver and solve system
A_petsc = as_backend_type(A)
b_petsc = as_backend_type(b)
x_petsc = as_backend_type(x.vector())
solver.set_operator(A_petsc)
solver.solve(x_petsc, b_petsc)
=============== *snap* ===============
What you might have missed is that the system needs to be consistent
(which is enforced in the above example, but may be true automatically
in other contexts, e.g., Navier-Stokes).
--Nico
On Wed, Jan 22, 2014 at 3:43 PM, Jan Blechta <[email protected]> wrote:
> On Wed, 22 Jan 2014 15:25:50 +0100
> Nico Schlömer <[email protected]> wrote:
>
>> Hi all,
>>
>> the problem is fully symmetric and can be solved with CG, even though
>> the null space is nontrivial (see, e.g., Iterative Krylov methods for
>> large linear systems, Henk A. van der Vorst). As a preconditioner you
>> can use AMG, the only thing that you have to take care of is the
>> coarse solve; by default, some LU strategy is deployed which messes
>> around with the null space. Take something iterative, e.g. Jacobi, and
>> you'll be fine.
>>
>> Check out
>>
>> prec = PETScPreconditioner('hypre_amg')
>> prec.parameters['hypre']['BoomerAMG']['relax_type_coarse'] = 'jacobi'
>> solver = PETScKrylovSolver('cg', prec)
>> solver.parameters['absolute_tolerance'] = 0.0
>> solver.parameters['relative_tolerance'] = tol
>> solver.parameters['maximum_iterations'] = 100
>> solver.parameters['monitor_convergence'] = verbose
>> # Create solver and solve system
>> A_petsc = as_backend_type(A)
>> b_petsc = as_backend_type(b)
>> phi_petsc = as_backend_type(phi.vector())
>> solver.set_operator(A_petsc)
>> solver.solve(phi_petsc, b_petsc)
>
> This does not seem to work. Have you tried it? Can you post a complete
> code? In my opinion, you can't solve it using CG because the problem is
> indefinite. At least, without some additional trick.
>
> Jan
>
>>
>> Alternatively, the PETSc option can be set using the PETScOption
>> framework that Garth added a while ago.
>>
>> I've encountered the same problem a while ago in the pressure update
>> for Chorin's method (Navier-Stokes solvers).
>>
>> Cheers,
>> Nico
>>
>> On Tue, Jan 21, 2014 at 8:15 PM, Jan Blechta
>> <[email protected]> wrote:
>> > On Tue, 21 Jan 2014 11:01:50 -0800
>> > Nikolaus Rath <[email protected]> wrote:
>> >
>> >> On 01/21/2014 09:43 AM, Simone Pezzuto wrote:
>> >> > 2014/1/21 Jan Blechta <[email protected]
>> >> > <mailto:[email protected]>>
>> >> >
>> >> > On Tue, 21 Jan 2014 17:18:54 +0000
>> >> > "Garth N. Wells" <[email protected] <mailto:[email protected]>>
>> >> > wrote:
>> >> >
>> >> > > On 2014-01-21 17:01, Nikolaus Rath wrote:
>> >> > > > Hello,
>> >> > > >
>> >> > > > I noticed that the neumann-poisson demo
>> >> > > >
>> >> >
>> >> > (http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/documented/neumann-poisson/python/documentation.html
>> >> >
>> >> > <https://urldefense.proofpoint.com/v1/url?u=http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/documented/neumann-poisson/python/documentation.html&k=Izx05CQZXsnLXkTIfmT7FQ%3D%3D%0A&r=1KW6QPJUrZMjRkn7m6Ouj0V90HWobEY8fXrhlFmuc%2Bc%3D%0A&m=v%2F3Ze6UBI9gR2l%2Fu82%2F4Y4fHujXMqW47hgHfmNXSjGw%3D%0A&s=e188e7efbe5d6d1deabbd7fb92be8e35f91642b13aacc34896849ef0c60c62e7>)
>> >> > > > fails when using a different solver, e.g. when replacing
>> >> > > >
>> >> > > > solve(a == L, w)
>> >> > > >
>> >> > > > with
>> >> > > >
>> >> > > > solve(a == L, w,
>> >> > > > solver_parameters = {'linear_solver': 'cg',
>> >> > > > 'preconditioner': 'ilu'})
>> >> > > >
>> >> > >
>> >> > > This problem needs very careful treatment when using
>> >> > > iterative solvers. Simple block-box preconditioners and
>> >> > > solvers will very likely fail.
>> >> >
>> >> > AMG preconditioning based on operator
>> >> >
>> >> > (inner(grad(u), grad(v)) + c*d)*dx
>> >> >
>> >> > could perform well. This operator does not have a dense row
>> >> > like the original one. This is a strategy similar to
>> >> > demo_stokes-iterative.
>> >> >
>> >> >
>> >> > In this case the preconditioner is singular (pure neumann), no it
>> >> > cannot be used.
>> >> >
>> >> > As Garth was mentioning, this problem is delicate for iterative
>> >> > solver, not only because
>> >> > its indefiniteness, but because the Lagrangian constraint you're
>> >> > imposing yields
>> >> > a column (the last one) of the full matrix that belongs to the
>> >> > kernel of the top-left block.
>> >> >
>> >> > Since the nullspace is at hands, I would provide it to the solver
>> >> > and then use CG+AMG,
>> >> > with Jacobi relaxation at coarser scale instead Gauss elimination
>> >> > (at least with petsc boomeramg).
>> >>
>> >>
>> >> Why is there a nullspace? Doesn't the \int u = 0 constraint remove
>> >> the remaining degree of freedom resulting from the pure Neumann
>> >> boundary conditions?
>> >
>> > It does not remove any DOF. It adds just one DOF - the Lagrange
>> > multiplier and an equation which makes the system regular. But this
>> > is evil as corresponding row and column in the assembled matrix is
>> > dense/full.
>> >
>> > Nullspace is an alternative (and probably better) approach.
>> >
>> > Jan
>> >
>> >>
>> >>
>> >>
>> >> Thanks,
>> >> Nikolaus
>> >>
>> >> _______________________________________________
>> >> fenics mailing list
>> >> [email protected]
>> >> http://fenicsproject.org/mailman/listinfo/fenics
>> >
>> > _______________________________________________
>> > fenics mailing list
>> > [email protected]
>> > http://fenicsproject.org/mailman/listinfo/fenics
>
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
