Diskontinuierliche Galerkin / Poisson / Fenics

Ich versuche, die 2D-Poisson-Gleichung mit der diskontinuierlichen Galerkin-Methode (DG) und der folgenden Diskretisierung zu lösen (ich habe eine PNG-Datei, darf sie aber leider nicht hochladen):

Gleichung:

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

Neue Gleichungen:

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

Schwache Form mit numerischen Flüssen T und 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$$

Numerische Flüsse (IP-Methode):

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

mit

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

Ich habe ein einfaches Fenics-Python-Skript geschrieben, um die Gleichung zu lösen. Die Lösung, die ich bekomme, ist nicht gut. Ich würde mich sehr freuen, wenn jemand, der mit der DG-Methode vertraut ist, einen kurzen Blick auf das folgende Skript werfen und mir sagen könnte, was ich falsch mache.

Die kontinuierliche Galerkin-Formulierung, die ich im Skript hinzugefügt habe, bietet eine gute Lösung.

Vielen Dank im Voraus.

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

Um Ihr Problem in FEniCS zu implementieren, müssen Sie die Integrale in Bezug auf Grenzen durch Integrale in Bezug auf Kanten ersetzen. Dies führt zu Sprüngen / Durchschnittswerten in den Testfunktionen, die Sie in Ihrer Implementierung völlig vermissen. Daher ist das System nicht invertierbar und Ihre Lösung sieht nicht richtig aus. Gleichung (3.3) in Arnold et. al. 2002 bietet Ihnen ein Werkzeug zum Umschreiben der schwachen 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$$

Hier ist die Vereinigung Ihrer Kanten und das gleiche ohne Grenzen.$\Gamma$$\Gamma^0$

Jetzt sind Ihre Flüsse einwertig, was bedeutet, dass Sie die Sprünge Ihrer Flüsse fallen lassen können. Daher ist

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

Dies führt uns zu folgender Änderung Ihres Codes:

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

Ich hatte noch nicht die Zeit, dies tatsächlich zu versuchen, also sei dir möglicher Vorzeichenfehler usw. bewusst. Aber ich hoffe, das hilft trotzdem.

Literatur: DN Arnold, F. Brezzi, B. Cockburn, LD Marini: Einheitliche Analyse diskontinuierlicher Galerkin-Methoden für elliptische Probleme SIAM J. Num. Anal, 39 (2002), 1749 & ndash; 1779

Ja, mir hat wirklich etwas gefehlt!

Es funktioniert jetzt gut.

Vielen Dank für Ihre Hilfe!