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.
Thanks a lot in advance.
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] - 0.5, 2) + pow(x[1] - 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(q,grad(v))*dx \
- 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()