Implementing DG in FEniCS

8 months ago by
I am trying to solve the following problem:
  $-\left(K\left(x\right)u'\left(x\right)\right)'=f\left(x\right)$(K(x)u'(x))'=ƒ (x)  
I am using the discontinuous Galerkin method and imposing the boundary conditions weakly. I am trying to practice with a case where:
$u\left(x\right)=\left(1-x\right)e^{-x^2}$u(x)=(1x)ex2 ,   
$u'\left(x\right)=\left(2x^2-2x-1\right)e^{-x^2}$u'(x)=(2x22x1)ex2  , and
$f\left(x\right)=u''\left(x\right)=\left(-4x^3+4x^2+6x-2\right)e^{-x^2}$ƒ (x)=u''(x)=(4x3+4x2+6x2)ex2 .

However, the results I am getting are off from the analytical solution:

Here is the relevant code:

import numpy as np
from fenics import *
import matplotlib as mpl
import matplotlib.pyplot as plt

parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['optimize'] = True

L = 1.0
start = 0.0
h = 0.1
num_elements = int(np.ceil((L-start)/h))
mesh = IntervalMesh(num_elements,start,L)
V = FunctionSpace(mesh,"DG",1)

v = TestFunction(V)
u = TrialFunction(V)

u_left = 1.0
u_right = 0.0
bcs = Expression("(u_left-u_right)/(L-start)*(x[0]-start)+u_right",degree=1,
                 u_left=u_left, u_right=u_right, start=start, L=L)

K = Constant(1.0)
u_true = Expression("(1-x[0])*exp(-x[0]*x[0])", degree=1)
eps = 1.0 # Penalizes jumps in the solution.
# Penalizes jumps in the solution and the test function.
sig0 = Constant(1.0)
J0 = sig0/h*jump(v)*jump(u)*dS
# Penalizes jumps in the solution's and the test function's derivative.
sig1 = Constant(0.0)
J1 = sig1/h*jump(v.dx(0))*jump(u.dx(0))*dS
# The final weak form
f = Expression("exp(-x[0]*x[0])*(-4*x[0]*x[0]*x[0] + 4*x[0]*x[0]+ 6*x[0] - 2)", degree=1)
lhs_F = K*u.dx(0)*v.dx(0)*dx - avg(K*u.dx(0))*jump(v)*dS + eps*avg(K*v.dx(0))*jump(u)*dS + J0 + J1
rhs_F = f*v*dx - eps*K*bcs*v.dx(0)*ds + sig0/h*v*ds

u = Function(V)
ffc_options = {"optimize": True, "quadrature_degree": 2}
problem = LinearVariationalProblem(lhs_F,rhs_F,u,[],form_compiler_parameters=ffc_options)
solver  = LinearVariationalSolver(problem)

u_true_space = interpolate(u_true, V)
coor = mesh.coordinates()
x = coor[:,0]
u_calc = u.compute_vertex_values(mesh)
u_true_arr = u_true_space.compute_vertex_values(mesh)

Is this the correct way to implement a DG method in FEniCS?

Community: FEniCS Project
So many things can go wrong with human error in DG. Did you forget your boundary term in + sig0/h*v*ds

refer to the undocumented dg-poisson demo if you get stuck.

If, ultimately, you want to solve a nonlinear problem, checkout my side project.
written 8 months ago by Nate  
Hello Nate,

  Thank you for your reply and for mentioning your side project.

  I am trying to implement the following weak form in FEniCS. For now, you can assume all the h values are equal. Note that this is inspired by (but not identical to) the weak form in Chapter 1 of the this monograph (click here to see it):

\sum_{n=0}^{N-1} \int_{x_n}^{x_{n+1}} K(x)u'(x)v'(x)dx - \sum_{n=0}^{N}\Big\{ K(x_n)u'(x_n) \Big\}\Big[v(x_n) \Big] + \epsilon \sum_{n=0}^N \Big\{K(x_n)v'(x_n)\Big\} \Big[u(x_n)\Big] \\
+ \sum_{n=0}^N \frac{\sigma^0}{h_{n-1,n}} \Big[v(x_n)\Big]\Big[u(x_n)\Big] + \sum_{n=1}^{N-1} \frac{\sigma^1}{h_{n-1,n}} \Big[v'(x_n)\Big]\Big[u'(x_n)\Big] \\
= \int_{\Omega} f(x)v(x)dx - \epsilon\Big(K(x_0)v'(x_0)u(x_0) -K(x_N)v'(x_N)u(x_N)\Big) + \frac{\sigma^0}{h_{0,1}}v(x_0)u(x_0) + \frac{\sigma^0}{h_{N,N-1}}v(x_N)u(x_N)

Is my translation of this DG weak form to FEniCS code correct? If not, what would be the correct way?

Thank you for helping.
written 8 months ago by Saro  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »