Diskontinuierliche Galerkin / Poisson / Fenics

10

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:

(κT)+f=0

Neue Gleichungen:

q=κTq=f

Schwache Form mit numerischen Flüssen T und q :T^q^

qwdV=T(κw)dV+κT^nwdSqvdV=vfdV+q^nvdS

Numerische Flüsse (IP-Methode):

q^={T}C11[T]T^={T}

mit

{T}=0.5(T++T)[T]=T+n++Tn

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

Antworten:

10

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:

KThKqKnKϕKds=Γ[q]{ϕ}ds+Γ0{q}[ϕ]ds

Hier ist die Vereinigung Ihrer Kanten und das gleiche ohne Grenzen.ΓΓ0

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

KThKq^nKvKds=Γ0q^[v]ds+Ωq^nvdsKThKwnKκT^ds=Γ[w]κ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

Christian Waluga
quelle
Ja, mir hat wirklich etwas gefehlt.
Micdup
-2

Ja, mir hat wirklich etwas gefehlt!

Es funktioniert jetzt gut.

Vielen Dank für Ihre Hilfe!

micdup
quelle
2
Können Sie der Vollständigkeit halber beschreiben, was Sie vermisst haben und wie Sie es behoben haben?
Paul