### Implementing DG in FEniCS

494

views

0

I am trying to solve the following problem:

$-\left(K\left(x\right)u'\left(x\right)\right)'=f\left(x\right)$−(

$u\left(0\right)=1$

$u\left(1\right)=1$

I am using the discontinuous Galerkin method and imposing the boundary conditions weakly. I am trying to practice with a case where:

$K\left(x\right)=1.0$

$u\left(x\right)=\left(1-x\right)e^{-x^2}$

$u'\left(x\right)=\left(2x^2-2x-1\right)e^{-x^2}$

$f\left(x\right)=u''\left(x\right)=\left(-4x^3+4x^2+6x-2\right)e^{-x^2}$

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

$-\left(K\left(x\right)u'\left(x\right)\right)'=f\left(x\right)$−(

`K`(`x`)`u`'(`x`))'=`ƒ`(`x`)$u\left(0\right)=1$

`u`(0)=1$u\left(1\right)=1$

`u`(1)=1I am using the discontinuous Galerkin method and imposing the boundary conditions weakly. I am trying to practice with a case where:

$K\left(x\right)=1.0$

`K`(`x`)=1.0 ,$u\left(x\right)=\left(1-x\right)e^{-x^2}$

`u`(`x`)=(1−`x`)`e`^{−x2},$u'\left(x\right)=\left(2x^2-2x-1\right)e^{-x^2}$

`u`'(`x`)=(2`x`^{2}−2`x`−1)`e`^{−x2}, and$f\left(x\right)=u''\left(x\right)=\left(-4x^3+4x^2+6x-2\right)e^{-x^2}$

`ƒ`(`x`)=`u`''(`x`)=(−4`x`^{3}+4`x`^{2}+6`x`−2)`e`^{−x2}.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)
solver.solve()
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)
plt.figure()
plt.plot(x,u_true_arr,'b',label='Analytical')
plt.plot(x,u_calc,'ro',label='FEniCS')
plt.legend(loc=1)
```

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

Community: FEniCS Project

Please login to add an answer/comment or follow this question.

`+ 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.

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

\begin{multline}

\\

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

\end{multline}

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

Thank you for helping.