Galerkin discontinuo / Poisson / Fenics

Estoy tratando de resolver la ecuación de Poisson en 2D usando el método Discontinuous Galerkin (DG) y la siguiente discretización (tengo un archivo png pero no puedo cargarlo, lo siento):

Ecuación:

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

Nuevas ecuaciones:

$$q = \kappa \nabla T\\\nabla \cdot q = -f$$

Forma débil con los flujos numéricos T y Q :$\hat{T}$$\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$$

Flujos numéricos (método

$$\hat{q} = \{\nabla T\} – C_{11} [T]\\ \hat{T} = \{T\}$$

con

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

Escribí un sencillo script de Python para resolver la ecuación. La solución que obtengo no es buena. Realmente agradecería si alguien familiarizado con el método DG pudiera echar un vistazo rápido al script a continuación y decirme qué estoy haciendo mal.

La formulación continua de galerkin que agregué en el script brinda una buena solución.

Muchas gracias por adelantado.

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()

$$\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$$

$\Gamma$$\Gamma^0$

$$\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$$

Esto nos lleva a la siguiente modificación de su código:

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(q,grad(v))*dx \
  + dot( qhat, jump(v,n))*dS \
  + dot( qhatbnd, v*n)*ds

Todavía no tenía tiempo para intentarlo, así que tenga en cuenta los posibles errores de signos, etc. Pero espero que esto ayude de todos modos.

Referencias: DN Arnold, F. Brezzi, B. Cockburn, LD Marini: Análisis unificado de los métodos discontinuos de Galerkin para problemas elípticos SIAM J. Num. Anal, 39 (2002), 1749-1779

Sí, realmente me faltaba algo!

Está funcionando bien ahora.

¡Muchas gracias por tu ayuda!