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? Thanks, Nikolaus _______________________________________________ fenics mailing list [email protected] http://fenicsproject.org/mailman/listinfo/fenics
