Galerkin / Poisson / Fenics descontínuo

10

Estou tentando resolver a equação de Poisson 2D usando o método Descontinuous Galerkin (DG) e a seguinte discretização (eu tenho um arquivo png, mas não tenho permissão para carregá-lo, desculpe):

Equação:

(κT)+f=0

Novas equações:

q=κTq=f

Forma fraca com fluxos numéricos e \ hat {q} :T^q^

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

Fluxos numéricos (método IP):

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

com

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

Eu escrevi um script python fenics simples para resolver a equação. A solução que recebo não é boa. Eu realmente apreciaria se alguém familiarizado com o método DG pudesse dar uma rápida olhada no script abaixo e me dizer o que estou fazendo de errado.

A formulação contínua de galerkin que adicionei no script fornece uma boa solução.

Muito obrigado antecipadamente.

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
fonte

Respostas:

10

Para implementar seu problema no FEniCS, você deve substituir as integrais em termos de limites por integrais em termos de arestas. Isso introduz saltos / médias nas funções de teste, que você perde totalmente em sua implementação. Portanto, o sistema não é invertível e sua solução não parece correta. A equação (3.3) em Arnold et. al. 2002 fornece uma ferramenta para reescrever a forma fraca:

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

Aqui é a união de suas arestas e o mesmo sem limites.ΓΓ0

Agora seus fluxos são de valor único, o que significa que você pode soltar os saltos de seus fluxos. Portanto,

KThKq^nKvKds=Γ0q^[v]ds+Ωq^nvdsKThKwnKκT^ds=Γ[w]κT^ds

Isso nos leva à seguinte modificação do seu 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

Ainda não tive tempo para realmente tentar isso, portanto, esteja ciente dos possíveis erros de sinal etc. Mas espero que isso ajude de qualquer maneira.

Referências: DN Arnold, F. Brezzi, B. Cockburn, LD Marini: Análise unificada de métodos descontínuos de Galerkin para problemas elípticos SIAM J. Num. Anal, 39 (2002), 1749-1779

Christian Waluga
fonte
Sim, estava realmente faltando alguma coisa.
Micdup 28/05
-2

Sim, estava realmente faltando alguma coisa!

Está funcionando bem agora.

Muito obrigado pela sua ajuda!

micdup
fonte
2
Por uma questão de integridade, você poderia descrever o que estava faltando e como o corrigiu.
Paul