Iterative method to solve Poisson problem with discontinous Galerkin formulation

13 months ago by
Hi everyone,

I would like to solve the Poisson problem (with piecewise constant coefficient, homogenous Dirichlet boundary conditions and a constant source term), in a discontinous Galerkin formulation and using an iterative method. I have found that using a direct solver, FEniCS is not able to solve it if the mesh exceeds 2562 nodes, that is why I need an iterative method.

Usually, to solve the Poisson equation in its classical weak formulation, it is enough to use a gmres method with an ilu preconditionner, but it does not seem to work with the discontinous Galerkin formulation.

Then, my question is the following : which are the best iterative method and preconditionner to solve the Poisson equation in a discontinous Galerkin formulation ?

For now, the best solution I have found, testing every solver and preconditionner I know, was to use a Krylov solver with a minres method and no preconditionner. However, the solver converges slowly and is not able to reach a tolerance below 10-4 .

Here is a minimum working code :

N = 512
mesh = UnitSquareMesh(N, N)
V2 = FunctionSpace(mesh, 'DG', 2)

# Define class marking Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
# Define test and trial functions
u = TrialFunction(V2)
v = TestFunction(V2)
w = Function(V2)
# Define normal vector and mesh size
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
# Define the source term f, and coeffecient mu
f = Constant(1.)
mu = Expression("1.+(x[0]>=0.5)",degree=0)
# Mark facets of the mesh
boundaries = FacetFunction('size_t', mesh, 0)
DirichletBoundary().mark(boundaries, 1)
# Define outer surface measure aware of Dirichlet and Neumann boundParies
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Define parameters
alpha = 4.0
gamma = 8.0
# Define variational problem
a = mu*dot(grad(v), grad(u))*dx \
    - dot(avg(mu*grad(v)), jump(u, n))*dS \
    - dot(jump(v, n), avg(mu*grad(u)))*dS \
    + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
    - dot(mu*grad(v), u*n)*ds(1) \
    - dot(v*n, mu*grad(u))*ds(1) \
    + (gamma/h)*v*u*ds(1)
L = v*f*dx
# Compute solution
if N>256:
    # Assemble system
    A, bb = assemble_system(a, L)
    # Assemble preconditioner system
    b = inner(grad(u),grad(v))*dx
    P, btmp = assemble_system(b, L)
    # Create Krylov solver
    solver = KrylovSolver("minres", "none")
    solver.parameters['absolute_tolerance'] = 1.E-4
    solver.parameters['maximum_iterations'] = 100000
    solver.parameters['error_on_nonconvergence'] = True
    # Associate operator (A) and preconditioner matrix (P)
    solver.set_operators(A, P)
    # Solve
    solver.solve(w.vector(), bb)
    solve(a == L, w)

Thank you for your help.


Community: FEniCS Project
I should precise that I seek the solution in a finite element space of order 2. May be that is why the convergence of the solver is so slow...
written 13 months ago by Fabien Vergnet  

3 Answers

13 months ago by
Discontinuous Galerkin methods can be difficult to precondition. For a continuous method for Poisson, GMRES + ILU is not a good solver, The 'best' is multigrid preconditioning.  You would need to do some research to find out what the best preconditioner for the DG case.

If your use the PETSc parameter system, e.g. 
PETScOption.set("option_foo", "bar")​

or petsc4py you will have the full range of PETSc solvers and preconditioners at your disposal.

Thank you for your response Garth. I will keep looking for a solution and post it here if I find something usefull.

written 13 months ago by Fabien Vergnet  
13 months ago by
Chris Richardson's weak scaling demo has a decent PETSc setup for Poisson solved with AMG preconditioning and CG elements (not DG) that might be a good starting place:

As Garth says though, that may not be sufficiently good for DG. Post back and let us know how you get on.
13 months ago by
This won't be your ideal answer. I just added a linear Poisson demo to the the dolfin_dg library. It uses a conjugate gradient iterative solver and hypre AMG to precondition the linear system. Perhaps you can use it to figure out what's going on in your own formulation.

dolfin_dg / demo /
Hi Nate. Thank you for your response.
I tried to run your demo in the dolfin container, but I got the error "no module name dolfin_dg".

Do you know what the problem is?

written 13 months ago by Fabien Vergnet  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »