On 01/21/2014 11:15 AM, Jan Blechta 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.
I'm probably just using different terminology, but how can the system be regular if it has a nullspace? If there is u such that A.u = 0, I would say that A is singular, not regular. Best, Nikolaus _______________________________________________ fenics mailing list [email protected] http://fenicsproject.org/mailman/listinfo/fenics
