On Wed, 22 Jan 2014 15:56:51 +0100
Nico Schlömer <[email protected]> wrote:
> 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).
No, I missed that you was talking about a problem without Lagrange
multiplier, which is only semi-definite, not indefinite.
Regarding problem with multiplier, I tweaked your suggestion using
MINRES instead of CG and preconditioner
ap = (inner(grad(u), grad(v)) + c*d)*dx
Code is below. This converges much faster than your code. In fact, your
code can't converge to rel tol 1e-15, probably because 'the
consistency' is not met to such a precision.
To be honest - comparing the residuals of these two problems (w and w/o
multiplier) is not much meaningful. One should rather compare errornorms
using exact solution.
========================================
from dolfin import *
# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R
# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
ap = (inner(grad(u), grad(v)) + c*d)*dx
L = f*v*dx + g*v*ds
# Assemble
A = assemble(a)
Ap = assemble(ap)
b = assemble(L)
# Create Krylov solver
prec = PETScPreconditioner('hypre_amg')
prec.parameters['hypre']['BoomerAMG']['relax_type_coarse'] = 'jacobi'
solver = PETScKrylovSolver('minres', prec)
solver.parameters['absolute_tolerance'] = 0.0
solver.parameters['relative_tolerance'] = 1e-15
solver.parameters['maximum_iterations'] = 100
solver.parameters['monitor_convergence'] = True
# Solve system
w = Function(W)
solver.set_operators(A, Ap)
solver.solve(w.vector(), b)
(u, c) = w.split()
# Plot solution
plot(u, interactive=True)
========================================
Jan
>
> --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
