Galerkin discontinu / Poisson / Fenics

J'essaie de résoudre l'équation de Poisson 2D en utilisant la méthode de Galerkin discontinu (DG) et la discrétisation suivante (j'ai un fichier png mais je ne suis pas autorisé à le télécharger, désolé):

Équation:

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

Nouvelles équations:

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

Forme faible avec flux numériques et : $\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$$

Flux numériques (méthode IP):

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

avec

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

J'ai écrit un simple script python fenics pour résoudre l'équation. La solution que j'obtiens n'est pas bonne. J'apprécierais vraiment que quelqu'un qui connaît la méthode DG puisse jeter un coup d'œil au script ci-dessous et me dire ce que je fais mal.

La formulation de galerkin continue que j'ai ajoutée dans le script donne une bonne solution.

Merci beaucoup d'avance.

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

Pour implémenter votre problème dans FEniCS, vous devez remplacer les intégrales en termes de limites par des intégrales en termes d'arêtes. Cela introduit des sauts / moyennes dans les fonctions de test, ce qui vous manque complètement dans votre implémentation. Par conséquent, le système n'est pas inversible et votre solution ne semble pas correcte. Équation (3.3) dans Arnold et. Al. 2002 vous donne un outil pour réécrire la forme faible:

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

Ici, est l'union de vos bords et le même sans frontières.$\Gamma$$\Gamma^0$

Maintenant, vos flux sont à valeur unique, ce qui signifie que vous pouvez supprimer les sauts de vos flux. D'où

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

Cela nous amène à la modification suivante de votre code:

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

Je n'ai pas encore eu le temps d'essayer, alors soyez conscient des possibles erreurs de signe, etc. Mais j'espère que cela vous aidera quand même.

Références: DN Arnold, F. Brezzi, B. Cockburn, LD Marini: Analyse unifiée des méthodes de Galerkin discontinues pour les problèmes elliptiques SIAM J. Num. Anal, 39 (2002), 1749-1779

Oui, il me manquait vraiment quelque chose!

Cela fonctionne bien maintenant.

Merci beaucoup pour votre aide!