# Discontinuous Galerkin / Poisson / Fenics

I am trying to solve the 2D Poisson equation using the Discontinuous Galerkin method (DG) and the following discretization (I have a png file but I am not allowed to upload it, sorry):

Equation : $$\nabla \cdot( \kappa \nabla T) + f = 0$$

New equations : $$q = \kappa \nabla T\\\nabla \cdot q = -f$$

Weak form with numerical fluxes $\hat{T}$ and $\hat{q}$:

$$\int q \cdot w dV = - \int T \nabla \cdot (\kappa w) dV + \int \kappa \hat{T} n \cdot w dS\\ \int q \cdot \nabla v dV = \int v f dV + \int \hat{q} \cdot n v dS$$

Numerical fluxes (IP method): $$\hat{q} = \{\nabla T\} – C_{11} [T]\\ \hat{T} = \{T\}$$

with $$\{T\} = 0.5 (T^+ + T^-)\\ [T] = T^+ n^+ + T^- n^-$$

I wrote a simple fenics python script to solve the equation. The solution I get is not good. I would really appreciate if somebody familiar with the DG method could have have a quick look at the script below and tell me what I am doing wrong.

The continuous galerkin formulation that I added in the script gives a nice solution.

from dolfin import *

method = "DG" # CG / DG

# Create mesh and define function space
mesh = UnitSquare(32, 32)
V_q = VectorFunctionSpace(mesh, method, 2)
V_T = FunctionSpace (mesh, method, 1)
W = V_q * V_T

# Define test and trial functions
(q, T) = TrialFunctions(W)
(w, v) = TestFunctions(W)

# Define mehs quantities: normal component, mesh size
n = FacetNormal(mesh)

# define right-hand side
f = Expression("500.0*exp(-(pow(x - 0.5, 2) + pow(x - 0.5, 2)) / 0.02)")

# Define parameters
kappa = 1.0

# Define variational problem
if method == 'CG':
a = dot(q,w)*dx \
+ T*div(kappa*w)*dx \
+ div(q)*v*dx

elif method == 'DG':
#modele = "IP"
C11 = 1.

a = dot(q,w)*dx + T*div(kappa*w)*dx \
- kappa*avg(T)*dot(n('-'),w('-'))*dS \
\
- dot( avg(grad(T)) - C11 * jump(T,n) ,n('-'))*v('-')*dS

L = -v*f*dx

# Compute solution
qT = Function(W)
solve(a == L, qT)

# Project solution to piecewise linears
(q , T) = qT.split()

# Save solution to file
file = File("poisson.pvd")
file << T

# Plot solution
plot(T); plot(q)
interactive()


To implement your problem in FEniCS, you have to replace the integrals in terms of boundaries by integrals in terms of edges. This introduces jumps/averages in the test functions, which you entirely miss in your implementation. Hence, the system is not invertible and your solution does not look right. Equation (3.3) in Arnold et. al. 2002 gives you a tool to rewrite the weak form: $$\sum_{K\in\mathcal{T}_h}\int_{\partial K} q_K \cdot n_K \phi_K\,ds=\int_\Gamma [q] \cdot \{\phi\}\,ds + \int_{\Gamma^0} \{q\} \cdot [\phi]\,ds$$

Here $\Gamma$ is the union of your edges and $\Gamma^0$ the same without boundaries.

Now your fluxes are single-valued, which means that you can drop the jumps of your fluxes. Hence $$\sum_{K\in\mathcal{T}_h}\int_{\partial K} \hat{q}\cdot n_K v_K\,ds=\int_{\Gamma^0} \hat{q} \cdot [v]\,ds + \int_{\partial\Omega} \hat{q} \cdot n v\,ds\\ \sum_{K\in\mathcal{T}_h}\int_{\partial K} w\cdot n_K \kappa\hat{T}\,ds=\int_{\Gamma} [w] \cdot \kappa\hat{T}\,ds$$

C11 = 1.
qhat = avg(grad(T)) - C11 * kappa*jump(T,n)
qhatbnd = grad(T) - C11 * kappa*T*n

a = dot(q,w)*dx + T*div(kappa*w)*dx \
- kappa*avg(T)*jump(w,n)*dS \
- kappa*T*dot(w,n)*ds \
+ dot( qhat, jump(v,n))*dS \
+ dot( qhatbnd, v*n)*ds


I did not have the time yet to actually try this, so be aware of possible sign-errors etc. But I hope this helps anyway.

References: D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini: Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems SIAM J. Num. Anal, 39 (2002), 1749-1779

• Yes I was really missing something. – micdup May 28 '13 at 16:01

Yes I was really missing something!

It is working fine now.

Thank you very very much for your help!

• For completeness sake, could you describe what it was you were missing and how you fixed it. – Paul May 28 '13 at 16:41