On Tue, 21 Jan 2014 11:01:50 -0800 Nikolaus Rath <nr...@trialphaenergy.com> wrote:
> On 01/21/2014 09:43 AM, Simone Pezzuto wrote: > > 2014/1/21 Jan Blechta <blec...@karlin.mff.cuni.cz > > <mailto:blec...@karlin.mff.cuni.cz>> > > > > On Tue, 21 Jan 2014 17:18:54 +0000 > > "Garth N. Wells" <gn...@cam.ac.uk <mailto:gn...@cam.ac.uk>> > > 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 > fenics@fenicsproject.org > http://fenicsproject.org/mailman/listinfo/fenics _______________________________________________ fenics mailing list fenics@fenicsproject.org http://fenicsproject.org/mailman/listinfo/fenics
